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I  .   INTRODUCTION 

A .   GENERAL 

War  policies  and  plans  for  military  operations  made 
during   peacetime    are   significant   for   the   mission 
accomplishments  of  combat   operations   conducted  during 
wartime.   Since   experimentation   with   real   combat  is 
infeasible,   military   analysts   use   stochastic  combat 
simulation  models  to  study  the  effects  of  policy  making 
on  combat  operations.   The   analysts'   inferences  drawn 
from  the   results  of   these  models  are  important  to  the 
decision  maker   since  he   has  to   use  them   to  make  the 
same  decisions  about  military  operations  as  he  would  if 
he  could  experiment  with  real  combat  itself.  The  output 
data   from   these   models   are   realizations   of  random 
variables   distributed   around    the    values    of   the 
parameters    of    interest,    or    the    models'    true 
characteristics,   so   the   analysts   can   only  estimate 
these  parameters   with  error.   The  magnitude   of   error 
for  each  estimate  can  be  measured  in  terms  of  precision 
or   the   variance   of   the  estimate:  if  the  estimate  is 
unbiased  then  the  smaller  the  variance,  the  greater  the 
precision,   and    the   smaller   the   error.   Since   the 
decision  maker  wishes  to  make  decisions   that  are  based 
on   estimate(s)   with   a   quantifiable  error  bound,  the 


analysts  may  find  it  possible,  to  apply  specific 
statistical  techniques  to  measure  and  control  the 
variance  of  the  estimate(s)  to  obtain  a  prescribed  or 
at  least  quantifiable  level  of  precision. 

The  analyst's   capability  of  estimating  a  parameter 
of  interest  with  high   precision  depends   on  the  extent 
to   which   he   is   able  to  control  the  sample  variance. 
When  the  analyst  uses   the   mean   value   of   the  sample 
output  data   as  the  estimate  of  a  parameter  of  interest 
and   when   individual   samples    are   independent,   the 
coefficient  of  the  variance  of  the  estimator  is  reduced 
by  a  factor  of  1 /n  where   n  is   the  sample   size  of  the 
output   data.   A   large   sample  size  yields  an  estimate 
with  a  small  variance   and     high   precision.  Multiple 
replications  to   obtain  a   large  sample  size  in  complex 
stochastic  models   can   be   prohibitively   expensive  in 
terms  of   resources  like  money,  internal  computer  time, 
computer  storage  space,  etc.   This   is   especially  true 
for   large-scale,   complex  stochastic  combat  simulation 
models,  which  often  require   hours  rather   than  minutes 
for   a   single   computed   replication.   Since  available 
computer  time  is  a   compelling   constraint   on  military 
studies  competing   for  scarce  resources,  the  analyst  is 
usually  given  an  allocated   amount  of   time  to  simulate 
his  model.   This  specified   amount  of  time  may  affect  a 
desired  level  of  precision  of  the   estimate(s)  that  the 


analyst  wishes  to  obtain  from  the  simulation.  Since 
the  analyst  can  only  execute  a  fixed  number  of 
replications  within  this  block  of  time,  the  sample  size 
(number  of  replications)  may  not  be  large  enough  to 
achieve  a  variance  small  enough  to  give  the  analyst  an 
acceptable  precision  for  the  estimate(s).  Hence,  the 
analyst  must  either  accept  the  particular  level  of 
precision  and  error  associated  with  such  variance  or 
apply  other  specific  statistical  techniques  which  are 
more  likely  to  produce  a  smaller  variance,  and  hence  a 
level  of  precision  with  which  he  can  feel  more 
comfortable . 

An  economizing  scheme  in  simulation  to  reduce  the 
variance  of  the  estimator  is  to  intentionally  distort, 
control,  and  modify  the  random  properties  of  the  input 
variables  in  the  simulation  model.  The  output  data 
resulting  from  the  manipulation  of  these  random  numbers 
are  random  variables  which  are  designed  to  be  much 
closer  together  and  more  closely  distributed  around  the 
true  value  of  the  model's  parameter  of  interest  than  is 
the  case  with  simple  random  sampling.  A  sample 
distribution  resulting  from  such  a  variance-reducing 
scheme  has  the  same  mean  value  but  a  potentially 
smaller  variance  than  the  distribution  of  the  sample 
without  the  usage  of  this  scheme.  The  different 
techniques   for   doing   this  scheme  are  called  Variance 
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Reduction  Techniques  (VRTs).  The  effects  of  certain  of 
these,  when  applied  to  a  combat  model,  are  the  subject 
of  this  thesis. 

B.   BACKGROUND 

VRTs  were  initially  used  to  evaluate  multi- 
dimensional integrals.  They  have  since  been  applied  to 
small  Monte  Carlo  simulation  problems  but  have  not  been 
extensively  utilized  in  large  complex  stochastic 
simulation  models.  The  utilization  of  these  variance- 
reducing  techniques  in  real-world  combat  simulation 
models  is  even  less  common.  Consequently,  limited 
examples  of  applications  of  VRTs  in  these  simulation 
models  are  found  in  the  literature.  The  major  reason 
for  this  is  because  the  performance  of  the  VRTs  is 
suspected  to  be  uncertain  and  unpredictable.  The 
analyst  has  no  guarantee  that  the  usage  of  VRTs  will 
work  all  the  time.  Futhermore,  he  has  no  way  to  know 
beforehand  how  much  variance  reduction  he  will  get  from 
the  application  of  VRTs  whenever  they  are  effective. 
However,  VRTs,  in  our  opinion,  promise  to  be  powerful 
and  effective  tools  in  simulation  if  the  issues  of 
their  performance  in  specific  simulations  are 
understood.  In  this  section  we  will  describe  the  effect 
that  they   can  have  on  simulation  studies.  In  Chapter  V 


we   illustrate   their   effectiveness   for   a  particular 
combat  model . 

The  effectiveness  of  a  VRT  may  be  measured  by  the 
relative  efficiency  of  the  simulation  in  obtaining 
estimate(s)  with  the  utilization  of  this  scheme,  to 
the  efficiency  of  a  simulation  under  the  same 
conditions  without  the  VRT.  Efficiency  as  Handscomb 
(1969, p.  253)  defines  it  is 

Efficiency  =  1  /  (Variance  *  Work).        (1) 

Here  "Work"  generally  refers  to  computing  time. 
According  to  Handscomb,  variance  reduction  succeeds  if 
the  VRT  increases  efficiency.  From  Equation  1 ,  we  see 
that  a  decrease  in  variance  and/or  work  will  increase 
efficiency.  Hence,  variance  reduction  in  simulation  is 
more  than  solely  a  decrease  in  the  variance  of  the  mean 
of  the  estimators.  Handscomb( 1 969 ,p .  253)  calls  a 
technique  variance-reducing  if  it  "reduces  the 
variance  proportionately  more  than  it  increases  the 
work  involved"  or  "does  not  reduce  the  variance  at  all 
in  the  usual  sense,  provided  that  it  saves  work."  The 
work  involved  in  attaining  estimates  by  simulation  has 
many  attributes.  Hammersley  and  Handscomb( 1 964 ,p .  22) 
suggest  that  the  number  of  simulation  runs  epitomizes 
this  work.  However,  we   can   easily   measure   this  same 
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work  in  terms  of  computing  cost  or/and  simulation  time. 
For  it  is  the  availability  of  these  factors  that 
ultimately  determine  the  precision  of  the  estimators. 
Hence,  an  effective  VRT  may  not  only  produce  more 
precise  estimates  but  also  economize  the  time  and  costs 
associated  with  the  simulation  to  obtain  the  level  of 
precision  for  those  estimates.  The  efficiency  of  VRTs 
will  be  discussed  more  fully  in  Chapter  III. 

This  potential  saving  in  computer  time  has 
stimulated  the  interest  of  the  United  States  Army 
Concept  Analysis  Agency  ( CAA )  in  the  utilization  of 
VRTs.  CAA  has  studied  the  effectiveness  of  a  VRT  in  two 
of  its  larger  and  more  complex  simulation  models.  The 
results  of  these  studies  were  mixed  (Johnson,  Bates, 
and  Graham,  1985).  CAA  therefore  recommended  the 
continuation  of  the  studies  to  investigate  the 
applicability  of  VRTs  to  reduce  the  inherent 
variability  in  large,  complex,  stochastic  combat 
simulation  models. 

C.   PURPOSE  AND  OBJECTIVE 

The  purpose  of  this  thesis  is  to  provide  additional 
insight  into  the  applicability  of  VRTs  to  stochastic 
combat  models,  and  to  provide  a  base  for  future  studies 
in  the  application  of  VRTs  to  large-scale,  real  world, 
stochastic  combat  simulations.   The  objectives   of  this 
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thesis  are  plainly  to  identify  those  VRTs  that  are 
applicable,  and  then  to  exhibit  their  performance  in 
the  applications  to  a  class  of  stochastic  combat 
simulation  models.  The  question  to  be  answered  is:  "Can 
VRTs  be  identified  that  are  consistently  effective  for 
reducing  simulation  time  and  cost?" 

D.  PROBLEM 

The  problem  for  this  thesis  is  to  increase  the 
efficiency  of  a  stochastic  combat  simulation  model 
utilizing  VRTs  in  terms  of  (1)  increased  precision  of 
the  model's  estimates  for  an  allocated  amount  of 
simulation  time,  and  (2)  reduced  computer  time  for  a 
predetermined  level  of  precision. 

E.  APPROACH 

In  his  doctoral  dissertation,  Andreasson  (1972) 
showed  that  variance  reduction  in  queuing  systems  is 
influenced  by  (i)  the  transformation  of  random  numbers, 
(ii)  the  structure  and  parameters  of  the  simulation 
model,  and  (iii)  the  choice  of  the  model  response 
quantity.  Condition  (i)  is  an  attribute  of  the  VRTs. 
Conditions  (ii)  and  (iii)  are  characteristics  of  the 
model.  To  solve  the  problem  stated  above,  we 
investigate  the  effects  of  the  parameters  of  a 
stochastic  combat  model,   described   in   Chapter   IV  of 
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this  thesis,  on  variance  reduction.  We  then  use  those 
results  to  formulate  our  approach  to  increase  the 
efficiency  of  this  model  in  terms  stated  in  the  problem 
above . 

F.   ORGANIZATION  OF  THIS  THESIS 

This  thesis  is  organized  into  6  chapters.  Chapter  I 
is  the  Introduction  chapter.  Chapter  II  reviews  the 
literature  of  VRTs  in  simulation.  Chapter  III  discusses 
ways  of  measuring  the  efficiency  of  a  VRT  and  explores 
the  tradeoffs  of  measuring  for  increased  precision  of 
estimation  and  reduced  computer  time.  Chapter  IV  is 
concerned  with  the  simulation  of  a  stochastic  combat 
model  and  the  programming  for  variance  reduction  in  the 
simulation  model.  Chapter  V  deals  with  the 
applicability  and  performance  of  VRTs  in  the  simulation 
model.  In  Chapter  VI,  we  make  conclusions  about  the 
applicability  of  VRTs  in  stochastic  combat  models  and 
provide  recommendations  about  their  use  in  larger  and 
more  complex,  stochastic  models  that  are  used  to  study 
real-world  combat  systems. 
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II .   REVIEW  OF  LITERATURE 

A.  INTRODUCTION 

The  VRTs  that  we  use  to  solve  the  problem  stated 
in  Chapter  I  of  this  thesis  are  antithetic  variates  and 
stratified  sampling.  But  first  we  review  the  literature 
of  variance  reduction  in  simulation.  This  chapter 
concentrates  on  the  practical  applications  of  VRTs  in 
simulation  models.  We  present  a  brief  summary  of  works 
of  scholars  and  experts  on  this  subject.  We  then 
describe  the  basic  concepts  of  two  VRTs  that  we  feel 
are  applicable  to  large-scale,  complex,  stochastic 
combat  simulation  models.  It  is  these  two  VRTs  whose 
performances  we  later  exhibit  in  the  combat  model  in 
this  thesis. 

B.  SUMMARY  OF  PREVIOUS  WORKS 

In  the  last  15  years  interest  in  VRTs  in  simulation 
has  stimulated  much  activity  on  this  topic  in  the 
Operations  Research  community.  This  section  does  not 
comprehensively  review  all  works  that  have  been  written 
in  the  literature,  but  it  presents  a  brief  overview  of 
the  utilization  of  VRTs  in  simulation.  The  purpose  of 
this  section  is  to  summarize  some  of  the  studies  of  the 
general  applicability  of  VRTs  in  simulation. 
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Hammersley  and  Handscomb  (-1964)  reviewed  many  of 
the  simplest  ideas  of  variance  reduction  in  simple 
Monte  Carlo  problems  as  they  can  be  applied  in  the 
fields  of  Mathematics  and  Physics.  Their  most  easily 
understood  examples  and  outstanding  successes  were  the 
evaluations  of  integrals  and  applications  to  particle 
physics.  Handscomb  (1969,  p. 252-262)  later  suggested 
that  VRTs  be  adapted  to  simulation.  He  acknowledged 
difficulties  in  predicting  the  effectiveness  of  the 
techniques  in  particular  situations,  but  he  did 
propose,  in  practice,  "...  to  proceed  by  more  or  less 
inspired  trial  and  error,  learning  by  experience  which 
tools  serve  one  best  [or  which  techniques  are 
effective]."  He  also  stated  that  it  may  be  much  harder 
to  tell  how  much  variance  reduction  may  occur  in  large 
and  complicated  simulation  problems.  These  issues 
remain  major  concerns  for  one  using  VRTs  in  large, 
complex,  stochastic  simulation  models. 

Moy  (1969,  pp.  263-288)  adapted  several  VRTs  to 
simulation  and  investigated  their  applicability  to 
queuing  systems.  He  concluded  that  VRTs  were  indeed 
capable  of  working  in  the  simulation  of  queuing 
systems.  Kleijnen  (1974,  Ch .  Ill),  who  has  written 
probably  the  most  comprehensive  and  most  referenced 
documentation  on  the  subject  of  VRTs  in  simulation, 
discussed  the   relevant  differences  between  sampling  in 
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Monte  Carlo  problems  and  sampling  in  stochastic 
simulation  models.  He  showed  that  VRTs  may  be  adapted 
to  accommodate  these  differences.  Kleijnen  also 
presented  a  detailed  description  and  critical  appraisal 
of  six  techniques  so  devised  for  utilization  in  the 
simulation  of  large  complex  systems.  These  VRTs  are 
stratified  sampling,  importance  sampling,  selective 
samp_l_i_n£,  control  variates,  antithetic  variates,  and 
common  random  numbers.  These  six  sampling  techniques 
have  become  the  most  well-known  and  popular  VRTs  in  the 
literature . 

Other  less-known  VRTs,  however,  have  been  applied 
to  simulation.  McGraft  and  Irving  (1974)  survey  some  18 
different  techniques  for  implementation  in  large  scale 
simulation  problems.  McGraft  and  Irving  include  a 
comprehensive  listing  of  the  characteristics, 
advantages  and  disadvantages,  and  criteria  for 
applicability  to  large  simulation  models,  and 
demonstrate  the  effectiveness  of  several  of  these 
techniques  with  a  military  simulation  application. 

Many  other  articles  and  papers  have  been  written  on 
the  subject  of  VRTs.  There  are  too  many  of  them  to  list 
in  this  thesis,  but  the  survey  ranges  from  specific 
techniques  to  more  general  methods  in  simulation 
experimentation.  Some  of  the  most  recent  papers  written 
about  the   general  applicability   of  VRTs  in  simulation 
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are  Nelson  (1985)  and  Cheng  (1986).  Textbooks  that 
illustrate  the  applications  with  simple  but  excellent 
examples  of  variance  reduction  in  simulation  are  Gaver 
and  Thompson  (1973,  Ch .  12),  Fishman  (1978,  Ch .  3),  Law 
and  Kelton  (1982,  Ch .  11),  Morgan  (1984,  Ch .  7),  and 
Br at ley,  Fox,  and  Schrage  (1987,  Ch .  2). 

C.   DESCRIPTION  OF  VRTs  USED  IN  THIS  THESIS 

Moy  (1969,  p.  263-288)  experimented  with  antithetic 
variates  and  stratified  sampling  and  showed  that  they 
are  indeed  capable  of  significantly  decreasing 
variability  in  the  simulation  of  simple  queuing 
problems.  Likewise,  we  wish  to  achieve  similar  results 
when  we  apply  them  to  the  simulation  of  our  stochastic 
combat  model.  We  do  this  in  Chapter  V.  In  this  section, 
we  discuss  the  underlying  conditions  and  fundamental 
concepts  in  the  applicability  of  antithetic  variates 
and  stratified  sampling  in  simulation. 
1 .   Antithetic  Variates 

The  method  of  antithetic  variates  is  relatively 
well-known  in  the  literature  of  variance  reduction  in 
simulation  (Kleijnen,  1974).  It  is  one  of  the  most 
useful  VRTs  because  of  its  simplicity  and  general 
applicability.  When  the  method  of  antithetic  variates 
is  used,  the  sampling  process  is  modified  by  the 
manipulation   of   random   numbers.    A   simulation   run 
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produces   a   response   from   the.  original   sequence  of 

random   numbers   (r>;   then,   a   second   simulation  run 

produces   an   antithetic   response  from  the  sequence  of 

the  complementary  random  numbers  { 1 -r > .  The   average  of 

the   two   responses   is   an   observation   on  the  sample 

output  data  of   the   stochastic   simulation   model .  The 

mean  value  of  this  sample  is  estimated  as  the  parameter 

of  interest. 

The  variance  of  this  estimate  is  reduced  if  the 

responses   of   the   first   and   antithetic  runs  of  each 

replication   are   negatively   correlated.   Besides   the 

interchanging  of   the  random   numbers  in   each  run,  two 

other   conditions    must   occur    to   produce   negative 

correlation  between  the  runs.  First,  each  response  must 

be  a  monotonic  function  of  its  respective  random  number 

stream;  that   is,  large  values  in  each  stream  of  random 

numbers  should  have  an  opposite  effect   on  the  response 

than   the   small   values,   and   vice   versa.  The  second 

condition  is  that  the   responses  to   the  events   in  the 

first  run   must  be   synchronized  with   the  responses  to 

the   events   in   the   antithetic   run.  Synchronization, 

defined  by  Kleijnen  (1974,  p.  193),  occurs 

...if  the  i  '  th  random  number  r-^  generates  [in  the 
stream  of  the  first  run]  a  particular  event  (e.g., 
arrival  of  customer  j)  then  in  the  antithetic  run 
(1  -  r-^)  should  generate  the  same  event  (i.e.,  not 
the  arrival  of  customer  j  '  where  j '  ^  j  and  not  a 
service  time ) . 
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If  the  antithetic  variate-s  methodology,  coupled 
with  required  conditions,  is  designed  into  the 
simulation,  then  the  average  of  the  two  negatively 
correlated  responses  will  tend  to  produce  an  estimate 
with  a  high  degree  of  precision.  That  is,  if  by  chance 
(r)  yields  a  response  above  the  value  of  the  true 
parameter  of  interest,  then  { 1 -r  }  should  yield  a 
response  below  the  value  of  the  true  parameter .  When 
these  responses  are  averaged,  the  deviations  between 
the  responses  and  the  true  parameter  approximately 
offset  each  other  resulting  in  relatively  small  net 
variability  in  the  output  data.  This  idea  can  be  shown 
mathematically.  Let  X-|  be  the  response  of  the  first 
run;  X2,  the  response  of  the  antithetic  run;  and  Y,  the 
average  of  X-|  and  X2. 

Y  =  (Xt  +  X2)  /  2 

VARy  =  1/4  *  {  VARX1  +  VARX2  +  2  *  C0VX1>X2  > 
=  1/2  *  (  1  +  C0RRx-|)X2)  *  VARX1 

Clearly,  a  negative  C0RRX1  X2  reduces  VARy  If 
C0RRX1  X2  equals,  or  is  close  to,  -1 ,  then  the  VARy  is 
mathematically  zero  or  very  close  to  it.  Hence,  the 
antithetic  sampling  is  designed  in  simulation  models  so 
that  the  correlation  between  the  pair  of  responses  is 
as  close  to  -1  as  possible. 
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Monotonicity  and  synahronizat ion  must  be 
designed  into  a  simulation  program  for  a  particular 
model  .  Kleijnen  (1974),  Law  and  Kelton  (1982),  and 
Bratley,  Fox,  and  Schrage  (1987)  are  excellent 
references  that  discuss  ways  to  do  this.  We  discuss  our 
design  to  achieve  these  two  conditions  for  antithetic 
variates  in  our  model  in  Chapter  IV.  As  stated  before, 
the  method  of  antithetic  variates  is  simple  to 
implement  and  requires  little  to  no  extra  computer 
time.  Because  of  simplicity  of  this  VRT ,  examples  of 
its  applications  are  illustrated  in  nearly  every 
textbook  that  considers  the  subject  of  VRTs. 
2 .   Stratified  Sampling 

The  stratified  sampling  technique,  discussed  in 
this  section  and  applied  to  the  simulation  model  in 
chapter  V  of  this  thesis,  is  a  different  version  of  the 
stratified  sampling  that  Moy ,  Kleijnen  and  other 
experts  on  VRTs  have  adapted  to  simulation.  Handscomb  ( 
1969,  p.  261)  calls  this  particular  version  of 
stratified  sampling  another  form  of  antithetic 
sampling.  Andrsasson  (1972,  p.  6)  refers  to  it  as  an 
antithetic  transformation.  Gaver  and  Thompson  (1973, 
pp.  585-586)  name  it  stratification  extending  an 
antithetic  idea.  It  is  indeed  stratification  in  that 
the  sampling  process  is  modified  so  that  the  range  of 
random  numbers   is  divided  into  two  or  more  strata  from 
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which  the  simulation  runs  produce-  responses.  It  has  the 
antithetic  flavor  in  that  the  responses  in  all  strata 
are  averaged  together  to  get  an  observation  which  is 
part  of  the  sample  output  data.  This  technique  is  also 
similar  to  antithetic  variates  in  that  its  estimator  is 
an  average  of  correlated  responses  (Gaver  and  Thompson 
1973,  p.  586).  Likewise,  this  estimator  tends  to  have  a 
smaller  variance. 

In  our  review  of  this  technique,  we  saw  no 
necessary  conditions,  like  those  for  the  antithetic 
variates,  for  this  technique  to  be  successful  in 
simulation.  The  design  of  stratified  sampling  into  our 
simulation  model  in  Chapter  IV  is  similar  to  that  one 
given  in  Gaver  and  Thompson  (1973,  p . 586 ) . 

D .   SUMMARY 

An  abundant  amount  of  material  has  been  written  on 
the  subject  of  variance  reduction.  Techniques  used  to 
reduce  the  variance  in  Monte  Carlo  problems  have  been 
adjusted  to  do  the  same  in  simulation  models.  The 
applications  of  VRTs  in  simulation  have  been 
illustrated  in  queuing  systems  and  simple  textbook 
problems  but  successful  applications  to  larger,  more 
complex,  real-world  stochastic  simulation  models  have 
not  been  so  amply  reported.  There  is  no  guarantee  that 
VRTs  will  work  spectacularly  for  every  situation  in  the 
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simulation,  and  when  they  do  work  it  is  necessary  to 
estimate  the  magnitude  of  the  variance  reduction.  Pilot 
tests  are  encouraged  to  help  resolve  these  issues. 
Antithetic  variates  and  modified  versions  of  stratified 
sampling  are  two  of  the  more  simple  and  easily  employed 
VRTs  and  will  be  applied  to  a  stochastic  combat  model 
in  Chapter  V. 
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Ill .  EFFICIENCY  OF  VARIANCE  REDUCTION 

A.   INTRODUCTION 

In  the  last  chapter   we  reviewed   some  studies  that 

involved  VRTs.   In  this   chapter  we  discuss  the  problem 

of   measuring   the   efficiency   of   a    VRT .   Comparing 

variances  of   a  parameter  of  interest  obtained  from  the 

simulations   with   and   without    the   use    of   a   VRT 

respectively,  on   an  ordinal   scale,  may   reveal  if  the 

VRT  works,  but  it  provides  little  information  about  how 

well  the   VRT  works.  Clearly,  a  quantitative  measure  is 

more  desirable.  Therefore,  the  manner  or  scale  on  which 

the  efficiency   of  a   VRT  is  measured  should  provide  as 

much  information   as  possible   on  the   performance  of  a 

VRT.   In   particular,   it   should  provide  at  least  some 

base  to  answering  the  following  questions: 

(i)    "Does  the  VRT  work?" 

(ii)   "If  so,  how  great  is  the   variance  reduction  in 
terms  of   increased  precision  for  estimating 
the  parameter  of  interest?" 

(iii)  "How  great  is  the  variance  reduction  in   terms 
of  simulation   time  saved  for  estimating  the 
parameter  of  interest?" 

(iv)   "What  are  the  tradeoffs,  if  any,  between  the 

potential  increase  in  precision  and  the  economy 
of  simulation  time  when  applying  VRTs?" 

In  the  next  section  we  examine  two  methods  that  are 

usually    used    in    the   literature   to   measure   the 
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efficiency  of  a  VRT .  We  evaluate  them  in  terms  of  how 
well  they  answer  the  questions  above.  In  the  third 
section,  we  offer  a  third  alternative  which  is  a  hybrid 
of  the  two  previous  methods  for  measuring  the 
efficiency  of  a  VRT.  This  third  method,  we  think, 
answers  all  four  questions  above  and  is  used  to  measure 
the  efficiencies  of  the  antithetic  variates  and  the 
stratified  sampling  techniques  whose  performance  is 
exhibited  in  this  thesis.  In  the  fourth  section  of  this 
chapter  we  show  how  to  use  the  third  method  of 
measuring  the  efficiency  of  a  VRT  to  obtain  the 
tradeoffs  between  increased  precision  and  reduced 
simulation.  The  last  section  is  a  summary  of  this 
chapter . 

B.   ASSESSMENT  OF  VARIANCE  REDUCTION 

In  the  literature  the  efficiency  of  a  VRT  is 
usually  measured  by  (1)  a  decrease  in  the  variance 
(Method  #1  )  or  (2)  the  relative  efficiency  of  a 
simulation  to  obtain  an  estimate  using  a  VRT  to  the 
efficiency  of  the  simulation  using  no  VRT  (Method  #2). 
Henceforth,  we  refer  to  a  simulation  without  the  use  of 
a  VRT  as  crude  simulation. 

Method  #1  is  well  defined  in  Kleijnen  (1974,  pp. 
106-107).  Kleijnen  uses  this  method  by  defining  the 
efficiency  of  a  VRT  as  a  percentage  of  reduction  in  varianc 
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Method  #1  =  (Var0  -  Var }   )  /  Var0)  *  \Q>®%  (2) 

where  Var0  and  Var  -|  are  variances  obtained  in  the  same 
amount  of  simulation  time  for  crude  simulation  and 
simulation  applying  a  VRT  respectively.  The  measure  of 
efficiency  of  a  VRT  which  Kleijnen  introduces  may  be 
interpreted  as  that  portion  of  the  variance  which  is 
not  achieved  by  crude  simulation  but  is  obtainable  in 
the  same  amount  of  simulation  using  a  VRT.  The  sign  of 
this  portion  determines  whether  the  VRT  increases  or 
decreases  the  precision;  a  positive  sign  reveals  an 
increase  and  a  negative  sign,  a  decrease.  The  magnitude 
of  the  portion  indicates  how  much  of  the  precision  is 
increased  or  decreased  respectively.  With  this  method 
we  can  also  see  that  the  VRT  has  an  identical  effect  on 
reducing  simulation  time  for  a  prescribed  level  of 
precision  as  it  does  on  increasing  precision.  Method  #1 
provides  answers  to  three  of  the  questions  stated  in 
the  last  section,  but  it  does  not  resolve  the  question 
of  tradeoffs  for  increased  precision  and  reduced  time 
in  a  simulation  using  VRT. 

McGrath  and  Irving  (1974,  p.  295)  use  Method  #2  to 
measure  the  efficiency  of  a  VRT.  They  initially  used 
this   method,   shown   as   Equation   (3),   to  equate  the 
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relative  advantage  gained  in  simulation  time  by  using  a 
VRT. 

Method  #2  = 
Varg/Var-]  *  (Simulation  Timeg )/( Simulation  Time-])   (3) 

where  subscripts  0  and  1  represent  crude  simulation  and 
simulation  applying  a   VRT   respectively.   The  relative 
efficiency  that   McGrath  and  Irving  used  to  measure  the 
efficiency  of  a  VRT   results  in   a  factor   by  which  the 
efficiency  of  a  simulation  is  increased  or  decreased  by 
using  a   VRT.  If   the  value   of  this   factor  is  greater 
than  one,   then  the   VRT  works;  otherwise,  it  does  not. 
The  magnitude  of  this  reduction  is  the  actual   value  of 
the  factor.   For  example,  if  the  value  of  the  factor  is 
5,  then  the  simulation  applying  the   VRT  can   obtain  an 
estimate   in   1 /5   the   simulation  time  required  by  the 
crude  simulation  for  the   same  precision   level.  Method 
#2  may  be  viewed   either  as  the  reduction  in  simulation 
time  when  both   simulations   are   to   achieve   the  same 
variance,   or   as   an   increase   in  precision  when  both 
simulations  are  run  for  the  same   amount  of   time.  This 
method,  like   Method  #1 ,   answers  only   the  first  three 
questions   proposed   in   the   first   section    of   this 
chapter . 
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C.   THE  HYBRID  METHOD 

Methods  #1  and  #2  measure  increased  precision  at  a 
fixed  simulation  time,  or  a  reduced  simulation  time  at 
a  fixed  level  of  precision.  They  do  not,  on  the  other 
hand,  measure  increased  precision  at  a  level  of  reduced 
simulation  time,  or  vice  versa;  nor  do  they  provide  a 
means  to  explore  such  a  possibility.  In  Chapter  I  we 
emphasize  that  variance  reduction  may  increase 
precision  and  reduce  simulation  time.  The  efficiency  of 
a  VRT ,  in  our  opinion,  should  reflect  both  effects  so 
that  we  can  explore  the  tradeoff  of  any  combination  of 
precision  and  simulation  time.  Method  #3  offers  such 
possibility  and  answers  all  four  questions  in  the 
introduction  section  of  this  chapter.  It  is  a  mixture 
of  Methods  #1  and  #2.  Method  #3  has  Kleijnen's  idea  of 
reduction  in  variance  and  McGrath  and  Irving' s  use  of 
relative  efficiency.  We  define  the  efficiency  of  a  VRT 
as  a  relative  efficiency  (RE),  as  shown  in  Equation  4, 
and  later  define  it  in  terms  of  increased  precision 
(IP)  and  reduced  time  ( RT ) . 

Method  #3  = 
Efficiency-!  /  Efficiency0  (4) 

where  Efficiencyg   and  Efficiency-)  are  the  efficiencies 
of  the  crude  simulation   and  simulation   applying  a  VRT 
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respectively.  The  definition  of  the  efficiency  of  a 
simulation,  identified  by  Equation  1,  is  the  inverse  of 
the  product  of  the  sampling  variance  of  the  parameter 
estimate  and  the  work.  Henceforth,  we  equate  work  to 
simulation  time,  which  is  the  total  time  of  the 
simulation  model  to  obtain  a  parameter  of  interest  and 
a  specified  variance.  Such  time  may  be  computed  as  n 
replications  times  T  (average)  simulation  time  per  run. 
If  k  runs  are  in  a  replication,  then  simulation  time 
equals  the  product  of  kn  runs  and  T  (average) 
simulation  time  per  run.  When  these  variables  are 
substituted  in  Equation  1 ,  the  efficiency  equation 
becomes  Equation  5: 

EFFICIENCY  =  1  /  (  Var  *  k  *  n  *  T)        (5) 

If  we  are  to  increase  the  efficiency  of  a 
simulation  using  a  VRT ,  then  we  must  attempt  to 
decrease  one  of  the  parameters  in  Equation  2.  The 
variable  k  runs  per  replication  is  a  constant  of  the 
VRT.  Specifically,  the  antithetic  variates  constant  k 
is  two;  stratified  sampling  constant  k  is  three  in  our 
study  (it  can  be  more);  and  for  no  VRT,  the  constant 
value  of  k  is  one.  The  variable  T  is  model  dependent; 
that  is,  its  value  depends  on  the  input  parameters  of 
the  model.   Attempts  to   decrease  this   variable  may  be 
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futile;  futhermore  Bratley,  Fox,  and  Schrage  (1987,  p. 
48)  point  out  that  relatively  little  can  be  done  to 
decrease  it.  We  discussed  the  relationship  between  the 
variance  (Var)  and  replications  (n)  in  Chapter  I.  They 
are  , in  essence, the  variables  we  wish  to  decrease. 
Throughout  this  thesis  we  interchange  the  phrases 
decrease  in  variance  with  increased  precision  and 
reduction  in  replications  with  reduced  time 
(simulation).  If  we  substitute  Equation  5  into  Equation 
4,  we  get  Equation  6.  Note  since  T  (average)  simulation 
time  per  run  is  the  same  for  both  simulations,  it  is 
left  out  of  the  equation. 

RE  =  Varg/Var!  *  kg,/!^  *  n0/n1  (6) 

If  the  RE  value  in  this  equation  is  greater  than 
one,  then  the  VRT  successfully  increases  the  efficiency 
of  the  simulation  and  is  said  to  be  strong;  otherwise, 
it  is  said  to  be  weak.  A  strong  VRT  decreases  the 
variance  so  that  precision  is  increased  and  simulation 
time  is  reduced.  A  weak  VRT,  on  the  other  hand,  does 
not  decrease  the  variance  as  well  as  a  strong  VRT;  in 
fact,  a  very  weak  (or  subversive)  VRT  may  increase  the 
variance,  which  causes  a  reduction  in  precision  and 
necessitates  an  increase  in  simulation  time.  In  most 
simulation  models   a   VRT   may   be   strong   for  certain 
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conditions  and  weak  for  other,  conditions.  In  this 
thesis,  we  look  for  such  characterizations  of  the 
antithetic  variates  and  the  stratified  sampling  in  the 
stochastic  combat  model  we  describe  in  the  next 
chapter.  In  the  next  section,  we  define  Equation(6)  in 
terms  of  increased  precision  and  reduced  simulation 
time  . 

D.   TRADEOFFS  OF  GAINING  PRECISION  AND  SAVING  TIME 
Let  us   define  IP  as  a  decrease  in  variance, 

IP  -  (Var0  -  Vart)  /  Var0  (7) 

and  RT  as  a  reduction  in  simulation  time 


Then 


and 


RT  =-  (n0k0  -  n-^)  /  n0k0.  (8) 

Var-L  -  (1.0  -  IP)  *  Var0  (9) 

n-^  -  (1.0  -  RT)  x  n0k  (10) 


If  we  substitute  Equations  (9)  and  (10)  into 
Equation  (6),  we  get  Equation( 11 ) . 

RE  =  1/(1.0  -  IP)  -  1/(1.0  -  RT)         (11) 
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Equation  1 1  defines  the  relative  efficiency  which 
we  defines  as  Method  #3  of  measuring  the  efficiency  of 
a  VRT ,  in  terms  of  increased  precision  and  reduced 
time.  This  equation  resolved  the  unanswered  question 
identified  as  (iv)  in  the  introduction  section  of  this 
chapter.  With  this  equation,  we  can  examine  any 
combination  of  IP  and  RT .  For  example,  suppose  we 
measure  the  efficiency  of  a  VRT  to  have  a  RE  value  of  6 
for  the  same  amount  of  simulation  time  (Hint:  RT=0). 
Substituting  these  values  into  Equation  (11),  we  get  IP 
=  5/6  or  83  .3%  increased  precision. 

Suppose  we  only  need   to  increase   the  precision  to 
75%  instead  of  83.3%,  then  we  can  substitute  the  values 
for   lP=3/4   and   RE=6   (RE   should   not   change)   into 
Equation  11.   We  now  get  RT=1/3  (Note,  we   increase  the 
precision  75%  and  reduce   the   simulation   time  33.3%). 
Likewise,   with   RE=6   for   the   efficiency  of  the  VRT, 
examples   of   other   combinations   are  ( IP=2/3 , RT= 1 /2 ) ; 
( IP=1 /2,RT=2/3) ;  and  ( IP=0 ,RT=5/6 ) .  In  fact,  we  may  get 
any  combination  of  (IP,RT)   between   0   and   5/6.  Note, 
however,  if   we  want   to  increase   the  precision  beyond 
83.3%  or  5/6,  we  will   get   an   increase   in  simulation 
time.  That   is  the   tradeoff  in  terms  of  more  increased 
precision.  For   example,  we   will  have   to  increase  the 
simulation  time  to  2/3  or  66.7%  to  accommodate  an  IP  of 
90%  for  a  RE  value   of   6.   In   short,   the  information 
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obtained  from  method  #3  is  that  we  can  make  an  estimate 
more  precise  and  save  simulation  time  simultaneously. 

E .   SUMMARY 

In  the  literature,  there  are  generally   two  methods 
of  measuring   the  efficiency   of  a   VRT .  Method  #1  is  a 
decrease  in  the  variance;   Method   #2   is   the  relative 
efficiency    of    a    simulation   using   VRT   to   crude 
simulation.  Both  methods  may  determine  if   VRT  works  in 
a   simulation    model.   They    also   may   indicate   the 
magnitude  of  the  variance  reduction  in   terms  of  either 
increased   precision   for   a   fixed   simulation  time  or 
reduced  simulation  time  at  a  fixed   level  of  precision. 
In   this   chapter,   we   introduced   a   third   method  of 
measuring  the  efficiency   of   a   VRT.   It   is   a  hybrid 
between   Method   #1   and   Method   #2.   Method  #3  offers 
exploration  into  the  tradeoffs   of  increasing  precision 
and  saving   simulation  time  for  any  efficiency  value  of 
a  VRT. 
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IV.   SIMULATION  OF  A  STOCHASTIC  COMBAT  MODEL 

A.   INTRODUCTION 

In  the  last  chapter  we  discussed  how  to  measure  the 
efficiency  of  a  VRT  to  determine  how  much  we  may  save 
in  simulation  time  or/and  how  much  we  may  increase  the 
precision  of  a  parameter  obtained  by  crude  simulation. 
In  this  chapter  we  show  how  we  may  apply  VRTs  to  the 
simulation  of  a  combat  system.  The  combat  model  which 
we  have  chosen  to  simulate  and  to  apply  the  VRTs  is  the 
BCD  Markovian  model  developed  in  the  doctoral 
dissertation  of  Abdul-Latif  Rashid  Al-Zayani  (1986).  A 
modified  version  of  this  model,  formulated  by 
Professor  Donald  P.  Gaver ,  is  in  Appendix  A.  This 
stochastic  model  may  seem  very  simple,  but  its 
simulation  provides  invaluable  insights  into  the 
applicability  of  VRTs  to  stochastic  combat  model. 

Beside  being  stochastic,  the  BCD  Markovian  model  is 
also  discrete  and  dynamic  in  nature;  hence,  it  is  a 
discrete-event  simulation  model.  We  refer  those  readers 
who  want  to  know  about  the  nature  of  discrete-event 
simulations  to  simulation  textbook  such  as 
Morgan( 1 984 ) ,  Law  and  Kelton( 1 982  )  ,  or  Bratley,  Fox, 
and  Schrage( 1 987 ) .  In  this  thesis,  we  describe  the 
simulation   of   the   BCD   model   in   terms   of  discrete 
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events.  We  describe,  in  detail,  the  crude  simulation  of 
the  BCD  model  in  the  next  section.  In  Section  C,  we 
show  how  we  applied  antithetic  variates  and  stratified 
sampling  to  this  simulation.  We  summarize  the  chapter 
in  the  last  section. 

B.   CRUDE  SIMULATION  OF  THE  BCD  MARKOVIAN  MODEL 

We  discuss  the  crude  simulation  of  the  BCD  model  in 
four  parts.  First,  we  describe  the  combat  scenario; 
second,  we  define  the  characteristics  of  the  model; 
third,  we  explain  the  simulation  of  the  combat  process 
in  the  model;  and  finally,  we  discuss  a  FORTRAN 
simulation  program  written  for  the  model. 
1  .   The  Combat  Scenario 

As  part  of  an  air  defense  command,  a  wing  of 
aircraft  defenders  is  responsible  for  defending  an  area 
against  a  hostile  air  attack  from  a  group  of  bombers. 
When  detection  of  an  incoming  threat  occurs  a  flight  of 
D  defenders  is  launched  to  engage  B  bombers  making  the 
attack.  When  the  two  groups  are  within  aerial  combat 
range  the  defenders  seek  a  one-to-one  combat  engagement 
with  the  bombers  at  a  rate  i.  Only  one  free  defender 
can  engage  a  free  bomber  in  combat;  a  bomber  will 
generally  attempt  to  avoid  any  engagement  with  a 
defender.  A  combat  engagement  lasts  until  either  the 
bomber  is   killed  or  the  defender  is  killed.  A  defender 
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kills  a  bomber  at  a  rate  a,  and  a  bomber  kills  a 
defender  at  a  rate  9.  Hence,  at  any  instant  during  the 
combat  process,  a  defender  is  either  free  and  searching 
to  fight  a  bomber,  fighting  a  bomber,  or  killed.  A 
bomber  is,  likewise,  either  free  and  eschewing 
engagement  with  a  defender,  engaging  in  combat  with  a 
defender,  or  killed.  The  combat  process  is  continued 
until  either  force  is  completely  killed  off  or  the 
duration  of  the  combat  period  is  terminated  after  T 
units  of  time. 

2 .   The  Characteristics  of  the  BCD  Model 

Hartman  (1985,  p.  2-18)  characterizes  the 
structure  of  a  combat  simulation  model  as  combat 
entities ,  attributes ,  and  events .  We  use  these 
characteristics  to  simulate  the  combat  process  of  the 
BCD  model  in  the  next  subsection.  Combat  entities  in 
the  BCD  model  are  free  bombers,  free  defenders,  and 
combat  engagements.  Each  entity  has  attributes  that 
describe  a  combat  scenario.  For  the  bombers,  the 
attributes  are  the  number  of  bombers  and  the  rate  a 
bomber  shoots  down  a  defender;  for  the  defenders,  the 
number  of  defenders  and  the  rate  a  defender  shoots  down 
a  bomber;  and  for  the  combats,  the  number  of  combat 
engagements  and  the  rate  that  a  bomber  and  a  defender 
engage  in  combat. 
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Law  and  Kelton  (1982,  p . . 4 )  define  an   event  in 
a   discrete-event   simulation   as   an   occurrence  which 
changes  the  state  of  the  system.  The  BCD  model  has  five 
events.  The   first  event   is  the   initialization  of  the 
air   battle.   The   initialization    event   governs   the 
initial  battle   conditions.  The   next  three  events,  are 
the  interim  events  in  the  combat   process.  These  events 
are  (1 )  a  combat  between  a  bomber  and  a  defender,  (2)  a 
defender  killing  a  bomber,  and  (3)   a  bomber   killing  a 
defender.   The   occurrence   of  an  interim  event  changes 
the  state  of  the  combat  process  at  time  t.  The  state  of 
the  combat   process  of   the  BCD  model  is  represented  by 
the   tr ivar iate-Markov    process   { B( t ) , C( t ) , D( t ) ; t >0 > ; 
where,  B(t)   is  the   number  of   free  bombers  at  time  t, 
C(t)  is  the  number  of  combat  engagements  at  time  t,  and 
D(t)  is   the  number   of  free   defenders  at  time  t .  As  a 
Markov  process,  the  combat  process  moves  from   state  to 
state  according   to  one-step  transit  probabilities  that 
depend  only  on  the   current  state.   The  fifth   and  last 
event  in   the  combat   process  is  the  termination  of  the 
air  battle.  The  termination  event  manifests  the  "end  of 
the  battle"   conditions.  The  values  of  the  state  of  the 
system   at   the   occurrence   of   the   termination  event 
reflects   the   battle   outcome.   These   values   are  the 
numbers  of  bombers  and  defenders  that  are  alive   at  the 
end  of   this  air   battle.  We  will  consider  these  values 
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as  the  parameters  of  interest  in_ the  simulation   of  the 
BCD  model . 

3 .   Simulation  of  the  Combat  Process 

We  simulated  the  combat  process  by  maintaining 
a  "bookkeeping  account"  of  the  changes  in  the  state  of 
the  combat  process  as  the  events  occur.  The  process 
begins  in  the  initial  state  { B( t ) , C( t ) , D( t ) ; t =0 }  with 
the  initialization  event  being  the  commencement  of  the 
air  battle.  Henceforth,  we  let  a  value  of  B(t)  equal  b, 
a  value  of  C(t)  be  c,  and  of  D(t)  be  d.  The  interim 
events  change  the  value  of  the  state  of  the  combat 
process  as  following: 

EVENT  STATE 

New  combat  b-1 ,c+1 ,d-1 

Bomber  kills  Defender  b,c-1 ,d+1 

Defender  kills  Bomber  b  +  1  ,c-1  ,d 

The  combat  process  spends  X(b,c,d)  units  of 
sojourn  time  in  state  (b,c,d)  until  another  event 
occurs.  The  sojourn  time  X(b,c,d)  is  a  random  variable 
distributed  exponentially  with  mean 

P(b,c,d)  =  1/(iBD  +  (a  +  9)C) 
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for  the  time  to  the  next  event  in- the  combat  process  of 
the  air  battle,  given  that  at  time  t  the  state  was 
(b,c,d).  Equation  12  is  this  sojourn  time.  To  derive 
this  expression,  we  use  the  inverse  transform  method  to 
obtain  unit-mean  exponential  random  variables,  where  V-j 
is  the  j th  random  number  in  the  sequence  of  a  stream  of 
uniform  random  numbers.  The  inverse  transform  method  is 
discussed  in  the  simulation  textbooks  listed  in  the 
reference  section  of  this  thesis. 

X(b,c,d)  =  -  P(b,c,d)  *  ln(V-j)  (12) 

We  use  the  value  of  Equation  12  to  advance  the 

simulated   time   of   the   air   battle   as  indicated  by- 
Equation  1  3  . 

t  =  t  +  X(b,c,d)  (13) 

The  combat  process  moves  to  another  state  when 
another  event  occurs.  The  probability  of  a  specific 
interim  event  occurring  is  governed  by  an  embedded 
Markov  chain  whose  transition  probabilities  are  as 
follows : 
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EVENT  .   PROBABILITY 

New  Combat  *BD  *  P(b,c,d) 

Defender  kills  Bomber  aC  *  P(b,c,d) 

Bomber  kills  Defender  eC  *  P(b,c,d) 

We  again  use  the  inverse  transform  method  to 
obtain  the  conditions  for  an  interim  event  to  occur  and 
to  induce  the  change  in  the  state  of  the  combat 
process.  V-j  is  the  j  th  random  number  in  the  sequence  of 
a  different  stream  of  random  numbers.  These  conditions, 
events,  and  changes  in  the  state  of  the  combat  process 
are  listed  below. 

CONDITION  EVENT  STATE 

VjiiBD*P(b,c  ,d)     New  Combat  b-1,c  +  1,d+1 

V-j  >iBD*P(b  ,c  ,d)     Defender  kills  Bomber    b,c-1,d+1 

and 

Vj  <(iBD+aC)*P(b,c,d) 

otherwise  Bomber  kills  Defender    b+1 ,c-1 ,d 

Thus,  we  (i)  generate  a  uniform  random  number 
to  choose  which  interim  event  has  occurred,  (ii)  update 
the  state  of  the  combat  process,  (iii)  generate  another 
uniform  random  number  and  transform  it  to  an 
exponential  random  variable  to  determine  the  unit  of 
time  until  the  occurrence  of  the  next  event  and  to 
advance  the   simulated  time   of  the   combat  process.  We 
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repeat  this  procedure  until  the  occurrence  of  the 
termination  event.  The  termination  event  occurs  when 
(1)  all  the  bombers  are  killed  (B(t)  +  C(t)  =  0)  or  all 
the  defenders  are  killed  (D(t)  +  C(t)  =  0)  or  (2)  the 
time  duration  of  the  air  battle  has  expired  (t  >T 
units  of  time).  At  the  end  of  the  aerial  battle,  the 
combat  process  is  in  state  (b,c,d)  from  which  we  can 
compute  the  number  of  live  bombers  B  (B(t)  +  C(t))  and 
the  number  of  live  defenders  D  (D(t)  +  C(t)).  These  are 
values  of  random  variables  for  one  possible  battle 
outcome . 

4 .   The  FORTRAN  Simulation  Program 

We  coded  the  crude  simulation  of  the  BCD  model 
in  FORTRAN.  This  FORTRAN  program,  consisting  of  a  main 
program  and  four  subroutines,  is  in  Appendix  B.  The 
main  program  begins  in  an  interactive  mode.  The  program 
reads  the  values  for  the  attributes  of  a  combat 
scenario  from  the  terminal  and  sends  them  to  the  BATTLE 
subroutine.  BATTLE  runs  N  replications  of  the  combat 
process  and  returns  the  summary  statistics  of  the 
outcome  of  N  battles  to  the  main  program.  The  main 
program  sends  them  to  the  STAT  subroutine.  STAT 
analyzes  these  battle  statistics  in  terms  of  parameter 
estimates  and  returns  the  values  of  these  parameter 
estimates  to   the  main   program.  The   main  program  then 
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sends   these   values   of    parameter   estimates    to   a 
formatted  output  file. 

The  BATTLE  subroutine  calls  two  subroutines 
UNIFOR  and  EXPON  for  the  generation  of  U(0,1)  random 
numbers.  These  two  subroutines  implement  the 
congruential  pseudo-random  number  generator 

ui+1  =  168071^  mod  (2**31  -  1)  (14) 

discussed  and  tested  by  Lewis  and  Orav  (1985,  Ch .  V). 
UNIFOR  generates  a  sequence  of  uniform  random  numbers 
for  the  selection  of  the  occurrence  of  an  interim 
event.  EXPON  generates  a  sequence  of  uniform  random 
numbers  for  the  computation  of  the  unit  of  sojourn  time 
in  a  state. 

The  STAT  subroutine  performs  statistical  output 
analysis  for  the  simulation.  It  computes  the  means  and 
variances  of  the  sample  distributions  of  live  bombers 
and  defenders.  The  sample  means  for  bombers  and 
defenders 

B  =  SUM  Bt  /  N  (15) 

D  =  SUM  D-l  /  N  (16) 

are  unbiased   (point)  estimators  of  E[B(t)]  and  E[D(t)] 
respectively.  Similarly,  the  variances 
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VARg  =  SUM  (  B±    -  B)**2  /  N  *  (N-1)        (17) 
VARf  =  SUM  (  B±    -    D)**2  /  N  x  (N-1)        (18) 

are    unbiased    estimators     of     VAR[E[B(t)]]    and 
VAR[E[D(t)]]  respectively  (Larson,  1982). 

C.   PROGRAMMING  FOR  VARIANCE  REDUCTION 

In  Chapter  I,  we  noted  that  VRTs  modify  the 
sampling  of  random  numbers.  In  this  section,  we  discuss 
these  modifications  for  the  antithetic  variates  and 
stratified  sampling  in  the  simulation  of  the  combat 
process.  We  describe  the  changes  we  made  to  the  crude 
FORTRAN  simulation  program  for  the  simulations  using 
antithetic  variates  in  Section  1  and  stratified 
sampling  in  Section  2  respectively. 

1 .   Antithetic  Variates 

We  make  changes  to  the  subroutines  BATTLE, 
UNIFOR,  and  EXPON  of  the  crude  simulation  program  to 
use  the  antithetic  variates.  The  FORTRAN  program  for 
the  BCD  simulation  model  applying  antithetic  variates 
is  in  Appendix  C.  The  BATTLE  subroutine  computes  the 
values  of  the  parameters  for  one  replication  as  the 
average  values  of  the  battle  outcomes  from  a  pair  of 
runs  of  the  combat  process.  We  obtain  the  values  of  the 
battle  outcome  from  the  first  run  by  using   a  stream  of 
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uniform  random  numbers  (U)  and  those  of  the  second  run 
by  using  a  stream  of  complementary  random  numbers  { 1  - 
U).  Since  we  are  attempting  to  decrease  the  variance  of 
the  estimates  by  inducing  negative  correlation  between 
these  two  runs,  we  want  to  minimize  this  negative 
correlation.  We,  first,  induce  a  negative  correlation 
between  the  two  runs  by  creating  monotonicity  between 
the  random  numbers  and  the  values  of  the  battle  outcome 
within  each  run.  We  then  minimize  this  negative 
correlation  by  synchronizing  the  sequences  of  random 
numbers  {U>  and  the  complement  (1-U)  (Bratley,  Fox, 
Schrage  1987,  p. 47).  Kleijnen  (1974,  p . 1 87 )  shows  that 
a  random  variable  generated  by  the  inverse  transform 
approach  is  monotonic.  Hence  we  have  monotonicity  in 
the  simulation  since  we  used  the  inverse  transform 
method  to  generate  the  uniform  random  variables  in  the 
simulation  of  the  BCD  model. 

Law  and  Kelton  (1982,  p.  352)  indicate  that  the 
inverse  transform  approach  also  facilitates  the 
maintenance  of  synchronization  .  With  this  method,  we 
use  only  one  uniform  random  number  per  sequence  to 
obtain  the  desired  random  variable  for  each  event  in 
the  combat  process;  as  contrasted  with  other  methods, 
like  the  rejection  method,  where  we  may  use  many  random 
numbers  to  produce  a  single  value  for  the  desired 
random  variable   of  the   same  events.   Thus  we  initiate 
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synchronization  in  the  model  when  we  use  the  inverse 
transform  method;  nevertheless,  we  must  still  preserve 
it. 

We  risk  losing  this  synchronization  in   the  BCD 
model  because  the  number  of  interim  events  occurring  in 
the  combat  process  per  run  is  a  random   variable.  Hence 
the  number  of  events  occurring  in  the  combat  process  in 
the  first  run  may  not  be   the   same   as   the   number  of 
events    occurring    in    the   combat   process   in   the 
antithetic   run.   Consequently,   the   number   of  random 
numbers   needed  in  the  antithetic  run  generally  differs 
from  that   required  in   the  first   run.  This  phenomenon 
leads  to   the  random   number  (Uj)   in  the  first  run  not 
being  synchronized  with  the  random  number  {  1  -U-j  }  of  the 
antithetic   run     (Kleijnen   1974,   p.   193).  In  other 
words,   the  complement  of  the  j th  uniform  random  number 
{  1 -U-j  }   is   not   used   for   the  jth  event  in  the  combat 
process  in  the   antithetic   run.   We   are   not   able  to 
control   the   random   number   of   interim  events  in  the 
combat  process,  but  we   can   manage   the   way   in  which 
UNIFOR   and   EXPON   generate   uniform  random  numbers  so 
that  synchronization  is  maintained  between  the   pair  of 
runs  per  replication. 

We  used  the  suggestions  of  Law  and  Kelton 
(1982,  p. 352)  and  Bratley,  Fox,  and  Schrage  (1987, 
p. 47)  to   maintain  the  synchronization  that  the  inverse 
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transform  method  has  initiated  in  the  BCD  simulation 
model.  We  modify  subroutines  UNIFOR  AND  EXPON  to 
generate  separate  streams  of  random  numbers  (U)  and 
complements  { 1 -U >  before  simulating  a  pair  of  runs  of 
the  combat  process.  When  the  subroutine  BATTLE  calls 
UNIFOR  and  EXPON,  it  receives  from  each  a  two- 
dimensional  array  of  random  numbers,  where  the  first 
column  contains  the  stream  (U)  and  the  second  column 
contains  the  stream  (1-U).  Hence,  we  use  the  first 
column  for  the  first  run  and  the  second  column  for  the 
antithetic  run.  This  approach  guarantees  that  if  the 
jth  event  in  the  first  run  uses  (U-j),  then  the  j th 
event  in  the  antithetic  run  will  use  { 1 -U-s  } .  We  do 
waste  some  of  the  random  numbers  in  the  arrays,  but  we 
do  it  judiciously.  Since  the  number  of  random  numbers 
used  in  the  combat  process  is  a  random  variable,  we  use 
only  those  random  numbers  that  we  need  in  each  column 
and  throw  away  the  remaining  so  that  no  overlap  is 
possible  for  the  next  pair  of  runs.  As  a  result,  we 
maintain  synchronization. 

The  last  change  we  make  to  the  crude  simulation 
for  the  utilization  of  the  antithetic  variates  is  in 
the  subroutine  BATTLE.  The  subroutine  BATTLE  computes 
the  values  of  parameters  for  each  replication  as 
follows : 
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Bi  =  (B1i  +  B2i)  /  2  (19) 

where 


B'i  =  the  number  of  live  bombers  from  the  first  run 


1 
B2-^  =  the  number  of  live  bombers  from  the  second  run 


and 

D±  =  (D1i  +  D2i)  /  2  (20) 

where 
D1  j_  =  the  number  of  live  defenders  from   the  first  run 
D2^  =  the  number  of  live  defenders  from  the  second  run 
2 .   Stratified  Sampling 

As  we  stated  in  Chapter  II,  stratified  sampling 
resembles  the  antithetic  variates  procedures,  and  so  do 
the  changes  to  the  crude  simulation.  Hence  we  make 
changes  similar  to  those  in  the  simulation  using 
antithetic  variates  for  the  simulation  using  stratified 
sampling.  The  FORTRAN  program  for  the  BCD  simulation 
model  using  stratified  sampling  is  in  Appendix  D.  We 
modify  subroutines  UNIFOR  and  EXPON,  where  each 
generates  a  three-dimensional  array  of  uniform  random 
numbers  from  the  three  strata 

St  =  (0,1/3),      S2  =  (1/3,2/3),      S3  =  (2/3,1) 

before  simulating  three  runs  of  the  combat  process  per 
replication.  Note  this  does  not  need  to  be  limited  to 
3;  we  could  have  done  more. 
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Gaver  and  Thompson  (1973,  p . 586 )  describe  this 
approach.  First,  a  U(0,1)  random  number  u-|  -j  is 
generated  and  placed  in  the  first  row  and  first  column 
of  the  three-dimensional  array.  Next,  1/3  is  added  to 
the  value  of  u-|  -j  with  a  subtraction  of  one  if  needed  to 
get  u-|2  within  the  range  of  0  and  "•  •  U12  is  placed  in 
the  first  row  and  second  column  of  the  array.  Next,  1/3 
is  added  to  the  value  of  u-|2»  and  if  necessary 
subtracted  by  one,  to  get  u-13.  u-13  is  placed  in  the 
first  row  of  the  third  column  of  the  array.  If 
subroutine  BATTLE  calls  for  k  random  numbers,  then  k 
U(0,1)  random  numbers  are  generated,  and  the  procedure 
obtains  a  value  for  each  of  the  kx3  cells.  BATTLE  uses 
the  first  column  of  random  numbers  in  the  array  for  the 
first  run,  the  second  column  for  the  second  run,  and 
the  third  column  for  the  third  run  of  the  combat 
process . 

The  values  of  the  parameters  for  each 
replication  are  the  average  values  of  the  battle 
outcomes  from  the  three  runs.  The  subroutine  BATTLE 
computes  the  values  of  these  parameters  as  follows: 

B±    =  (B1i  +  B2i  +  b5±)  /  3  (21  ) 

where 

B1  j_  =  the  number  of  live  bombers  from  the  first   run 
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B2i  =   the  number  of  live  bombers  from  the  second  run 
B^  =  the  number  of  live  bombers  from  the  third  run 
and 

Dt  =  (D1i  +  D2i  +  D^)  /  3  (22) 

where 

D1  -j_  =  the  number  of  live  defenders  from  the  first  run 

D2-^  =  the  number  of  live  defenders  from  the  second  run 

D^  =  the  number  of  live  defenders  from  the  third  run 

D .   SUMMARY 

The  simulation  of  the  BCD  model  is  a  discrete-event 
simulation.  It  begins  with  the  initialization  event  and 
ends  with  termination  event.  The  simulation  of  the 
combat  process  involves  generating  a  sequence  of  U(0,1) 
random  numbers  to  select  interim  event  occurrences  with 
changes  in  the  state  of  the  process  and  generating 
another  sequence  of  U(0,1)  random  numbers  to  determine 
the  unit  of  time  until  the  next  event  occurs  and  to 
advance  the  simulated  time  of  the  combat  process. 

The  programming  of  the  antithetic  variates  and 
stratified  sampling  modifies  crude  simulation. 
Monotonicity  and  synchronization  are  required  in 
generating  the  uniform  numbers  for  the  simulation  using 
these  VRTs .  Generating  random  numbers  by  the  inverse 
transform  method  guarantees  monotonicity.  Generating 
sufficient   random   numbers   by   the   inverse  transform 
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method   and   in   multi-dimensional  arrays  initiates  and 
maintains  synchronization. 
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V.   APPLICATIONS  OF  THE  ANTITHETIC  VARIATES  AND 

STRATIFIED  SAMPLING 


A.   INTRODUCTION 

We  are  now  prepared  to  demonstrate  the  application 
of  the  variance  reducing  techniques  to  the  simulation 
of  a  combat  stochastic  model.  In  this  chapter  we 
illustrate  the  performance  of  the  antithetic  variates 
(AV)  and  stratified  sampling  ( SS )  in  the  simulation  of 
the  BCD  model .  In  Chapter  IV  we  stated  that  the  mean 
and  variance  of  the  parameters  of  interest  estimated 
from  simulation  are  used  to  analyze  the  output  data  of 
the  model.  Usually  the  estimated  mean  is  of  primary 
interest  to  decision  makers,  and  the  estimation  of  the 
variance  is  secondary.  Since  we  use  the  variance  of  the 
parameters  estimated  from  the  simulation  of  the  BCD 
model  to  exhibit  variance  reduction,  we  will, 
henceforth,  focus  on  the  variance. 

We  examine  the  applicability  of  AV  and  SS  in  the 
BCD  model  by  simulating  many  scenarios  of  the  air 
battle  and  recording  increases  in  simulation 
efficiency.  We  investigate  AV  and  SS  performance  by 
mapping  a  response  surface  that  characterizes  the 
efficiency  of  variance  reduction  in  the  model.  In  the 
next  section  we  specify   the  scenarios   and  discuss  the 
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application  of  using  AV  and  SS  in  the  simulation  of 
these  scenarios.  We  then  build  response  models  that 
describe  the  performance  of  the  AV  and  SS  for  these 
scenarios  and  discuss  their  experimental  results  in 
Section  C.  We  present  a  brief  summary  of  the  chapter  in 
the  final  section. 

B.   APPLICABILITY  OF  ANTITHETIC  VARIATES  AND  STRATIFIED 
SAMPLING  TO  THE  BCD  SIMULATION  MODEL 

In  Section  B  of   the  previous   chapter  we  described 

the  general   scenario  of  the  BCD  model.  In  this  section 

we  specify  various  combat   scenarios  to   observe  how  AV 

and   SS   are   applicable   to   the  simulation  of  the  BCD 

model.  Recall  that  seven  attributes  characterize   a  BCD 

scenario.  Since   a  change  of  the  values  of  one  of  these 

attribute  will  produce  a   different  scenario,   we  chose 

to  change  the  values  for  three  attributes.  We  simulated 

10,  30  and  50  defenders  against  10,  30,   and  50  bombers 

at   "end   of   battle"   times   of   25,   75,   and  125.  We 

maintain  the  a  rate  for  a  defender  killing   a  bomber  at 

.05,  the  e  rate  for  a  bomber  killing  a  defender  at  .05, 

and  the  i    rate  for  a  bomber  and   a  defender   entering  a 

combat   engagement   at   .005.   We  also  initialize  every 

simulation   with   zero   combat   engagements.   Thus,   we 

observe  the   responses  of  the  simulations  at  3   "end  of 

battle"    times    in    9    different    scenarios.   This 
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arrangement    comprises    a    set    of   27   independent 
simulations . 

We  run  this  arrangement  for  crude  simulation, 
simulation  using  AV ,  and  simulation  using  SS.  As  a 
result,  we  perform  a  total  of  81  different  simulation 
experiments.  In  order  to  make  a  fair  assessment  of  the 
applicability  of  AV  and  SS ,  we  examine  the  variance 
obtained  from  the  same  amount  of  simulation  or  the  same 
numbers  of  simulated  battle  runs  for  every  simulation. 
We  run  90  battles:  this  equates  to  90  replications  in 
crude  simulation,  45  replications  in  simulation  using 
AV ,  and  30  replications  in  simulation  using  SS .  Table 
E.1  of  Appendix  E  contains  the  statistical  output  for 
crude  simulation;  similarly,  data  in  Tables  E.2  and  E.4 
are  from  simulations  using  AV  and  SS  respectively.  We 
use  Equations  6  and  7  to  measure  the  efficiency  of  the 
variance  reduction  (RE)  and  the  increase  in  precision 
of  the  parameters  from  the  simulation  (IP)  and  place 
the  AV  results  in  Table  E.3  and  SS  results  in  Table 
E.5. 

The  values  in  Tables  E.3  and  E.5  show  that  AV  and 
SS  respectively  are  applicable  in  the  BCD  simulation 
model.  A  RE  value  greater  than  unity  indicates  that  the 
VRT  is  effective  in  increasing  simulation  efficiency.  A 
positive  IP  value  exhibits  their  effectiveness  to 
increase  precision  of  the  desired  parameter.  With  these 
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two   values,   we   may   also  find  the  tradeoff  of  saving 
simulation  time  or   gaining   precision.   We  acknowledge 
that  these   values  are   all  values   of,  or  realizations 
of,  random   variables,  but   we  believe   that  the  tables 
show    that    the    variance   reduction   adheres   to   a 
stochastic   pattern.   That   is,   the   random   variables 
obtained  under   certain  scenarios  will  tend  to  have  the 
same   relationship   to   the   random   variables  obtained 
under   other   scenarios.   For   example,  the  data  in  the 
tables  suggest  that  high  RE  and   IP  values    correspond 
to  the  scenarios  that  start  combat  with  same  numbers  of 
bombers  and   defenders.  The   RE  and   IP  values  obtained 
under  these   scenarios  appear   consistently  higher  than 
the   values   under   all   other   scenarios.   Hence,   the 
variance  reduction   measured  by   these  RE  and  IP  values 
are  stochastically  greater  than   the  variance  reduction 
obtained   from   any   other   scenario.   Since   such  even 
combats  (i.e.  equal  combat   power)  are   inherently  more 
variable  in   outcome,  the   fact  that  variance  reduction 
is  greatest  there  is   certainly   welcome.   In   the  next 
section    we    attempt    to   conduct   a   more   thorough 
investigation   of   these   phenomena   so   that    we   may 
understand  how   the  AV  and  SS  perform  in  the  simulation 
of  the  BCD  model . 
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C.   PERFORMANCE  OF  ANTITHETIC  VARIATES   AND  STRATIFIED 
SAMPLING  IN  THE  BCD  SIMULATION  MODEL 

In  the  previous  section  we  saw  that  AV  and  SS  are 
applicable  to  the  simulation  of  the  BCD  model.  In  this 
section  we  examine  the  variability  of  uncertainty  in 
the  model  and  then  evaluate  the  applicability  of  AV  and 
SS  to  reduce  this  uncertainty.  We  explore  the  changes 
in  the  AV  and  SS  performance  and  examine  the 
relationships  of  factors  that  affect  these  changes. 
Results  of  this  analysis  reveal  the  characterization  of 
the  AV  and  SS  performance  in  the  BCD  model. 

1 .   Experimental  Design 

We  use  the  data  we  generated  in  Appendix  E  to 
fit  response  surfaces  that  describe  the  uncertainty  in 
the  values  of  parameters  in  the  BCD  model  and 
characterize  the  performance  of  AV  and  SS  over  the 
prescribed  range  of  values  in  the  three  factors: 
initial  numbers  of  Bombers  and  Defenders  and  "end  of 
battle"  Time.  We  code  the  three  factors  as 

x-|  =  (Time  -  75)  /  50, 
x2  =  (Defender  -  30)  /  20, 
X3  =  (Bomber  -  30)  /  20. 

Each   factor   has   3   levels.   Thus,   we   may   use  a  3^ 
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factorial    design    to      fit      the      data      with      the    response 
surface    equation 

E(y)=    eg    +    e-|X-|     +    p2x2    +    03x3    +    e-|1x1**2 
+    P22x2**2    +    033X3**2    +    e-|2x1*x2 
+    a-l^x^x^    +    B23x2*x3  (20) 

where 

E(y)=  expected  response 
00=  intercept 

p^=  linear  coefficient  for  factor  i 
pii=  Quadratic  coefficient  for  factor  i 
0ji=  interaction   coefficient  for  the 
interaction  of  factors  i  and  j 
x^=  level  of  factor  i. 
We  seek  to  obtain  the  maximum   information  from 
every  observation;   therefore,  we  chose  a  3^  Fractional 
factorial  design.  This  design  is  the  cuboctahedron  plus 
three   center   points   (John,   1971).   The  three  center 
points   provide   an   unbiased   estimate   of   error   and 
repeated  observations   which  permit  us  to  test  for  Lack 
of  Fit  of  the  response  surface  equation  we  obtained.  We 
use   the   cuboctahedron   design   to   fit  data  for  three 
response   surfaces:   (1)   variability    of   uncertainty 
inherent  in   the  battle  outcomes,  (2)  efficiency  of  AV , 
and  (3)  efficiency  of  SS .   This  design,   with  its  three 
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center  points,  and  the  data  to  which  it  is  used  to  fit 
are  shown  in  Tables  5.1,  5.2,  and  5.3.  The  variances 
data  in  Table  5.1  is  the  variability  of  uncertainty 
inherent  in  the  battle  outcome.  We  obtain  this  data 
from  the  variance  of  the  estimate  tabulated  in  the 
crude  simulation  table  in  Appendix  E.  The  RE  values 
data  in  Table  5.2  is  the  efficiency  of  AV  for  the 
estimation  of  the  defender  and  bomber  parameters.  We 
took  this  data  from  the  RE  values  in  Table  E.3.  In 
Table  5.3  is  the  RE  values  data  for  the  efficiency  of 
SS .  This  data  is  obtained  from  the  RE  values  in  Table 
E.5. 


TABLE  5.1  DESIGN  AND  VARIANCES  FOR  A  3  x  3  EXPERIMENT 
ON  THE  VARIABILITY  OF  UNCERTAINTY  INHERENT  IN  THE 
BATTLE  OUTCOME 


DECODED 

LEVELS 

CODED 

LEVELS 

VARIANCES 

Time 

Defen 

ier  Bomber 

x1 

x2 

x3 

Defender 

Bomber 

25 

1  0 

30 

—  1 

-1 

0 

.0328 

.0541 

125 

10 

30 

-1 

0 

.0036 

.  1725 

25 

50 

30 

- 1 

1 

0 

.  1687 

.0826 

125 

50 

30 

1 

0 

.3339 

.0359 

25 

30 

10 

- 1 

0 

- 1 

.0458 

.0350 

125 

30 

10 

0 

-1 

.2247 

.0045 

25 

30 

50 

- 1 

0 

.0628 

.  1  1  04 

125 

30 

50 

0 

.0362 

.3479 

75 

10 

10 

0 

-1 

- 1 

.0387 

.0358 

75 

50 

10 

0 

1 

- 1 

.  1427 

.0063 

75 

10 

50 

0 

-1 

.0064 

.  1654 

75 

50 

50 

0 

1 

.2717 

.2865 

75 

30 

30 

0 

0 

0 

.1618 

.  1895 

75 

30 

30 

0 

0 

0 

.0925 

.0966 

75 

30 

30 

0 

0 

0 

.  1693 

.  1589 
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TABLE  5.2  DESIGN  AND  RE  VALUES  FOR  A3  x  3  EXPERIMENT 
ON  THE  EFFICIENCY  OF  THE  ANTITHETIC  VARIATES  IN  THE  BCD 
MODEL 


DECODED  LEVELS 
Time  Defender  Bomber 


CODED  LEVELS 


x1 


x2 


x3 


RE  VALUES 


Defender   Bomber 


25 

1  0 

125 

1  0 

25 

50 

125 

50 

25 

30 

125 

30 

25 

30 

125 

30 

75 

1  0 

75 

50 

75 

10 

75 

50 

75 

30 

75 

30 

75 

30 

30 
30 
30 
30 
10 
10 
50 
50 
10 
10 
50 
50 
30 
30 
30 


0 
0 
0 
0 
0 
0 
0 


-1 

-1 
1 
1 

0 
0 
0 
0 

-1 
1 

-1 
1 

0 
0 
0 


0 
0 
0 


2 
1 

6 
7 
3 
4 
2 
1 


1 
1  0 

6 
13 


5 
5 
0 
2 
9 
0 
8 
8 


8.  1 
3.5 


1 
1 

4 
5 


5.4 
3.9 
2  .2 

1  .4 
3.0 


1 

4 
9 


10.9 
6.0 


1  0.8 


7.0 
7.0 
9.2 


TABLE  5.3  DESIGN  AND  RE  VALUES  FOR  A3  x  3  EXPERIMENT 
ON  THE  EFFICIENCY  OF  STRATIFIED  SAMPLING  IN  THE  BCD 
MODEL 


DECODED  LEVELS 


Time  Defender  Bomber 


CODED  LEVELS 


RE  VALUES 


x1 

x2 

x3 

—  1 

-1 

0 

-1 

0 

- 1 

1 

0 

1 

0 

- 1 

0 

- 1 

0 

- 1 

-  1 

0 
0 

0 

-1 

- 1 

0 

1 

- 1 

0 

-1 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Defender   Bomber 


25 

10 

125 

1  0 

25 

50 

125 

50 

25 

30 

125 

30 

25 

30 

125 

30 

75 

1  0 

75 

50 

75 

1  0 

75 

50 

75 

30 

75 

30 

75 

30 

30 
30 
30 
30 
10 
10 
50 
50 
10 
10 
50 
50 
30 
30 
30 


3 
2 
2 
2 
2 
1 
1 
1 

2 
3 
3 
2 
2 


6 
9 
0 
2 
0 
6 
9 
7 
6 
6 
0 
0 
7 
9 
5 


2 
2 
1 

1 

2 

1 

1 

2 

1 

1 

1 

3 

2 

3 

1 


1 

2 

9 

2 

8 

7 

4 

6 

5 

2 

5 

3 

7 

0 

5 
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2 .   Experimental  Analysis 

To  perform  the  statistical  analysis,  we  use  the 
Response  Surface  Regression  (RSREG)  procedure  in  the 
Statistical  Analysis  System  ( SAS )  computer  software 
package  on  the  IBM  370  mainframe.  With  this  procedure 
we  were  able  to  obtain  a  second  order  response-surface 
equation  by  least-square  regression,  check  for  model 
adequacy,  test  for  lack-of-fit,  and  identify  critical 
surface  values  which  were  useful  in  helping  to  describe 
the  shape  of  the  surface. 

We  fitted  Equation  20  to  the  data  in  the 
respective  design  tables  in  the  previous  section  and 
obtained  multiple  response  surface  equations  and 
multiple  analysis  of  variance  (ANOVA)  tables  for 
corresponding  responses.  We  assess  the  adequacy  of  each 
equation  and  test  for  fit  from  its  corresponding  ANOVA 
table.  From  each  response  surface  equation,  we 
generated  additional  data  to  obtain  contour  plots.  We 
plotted  contours  of  variability  of  uncertainty  for  the 
initial  numbers  of  bombers  and  defenders  at  various 
Times.  Similarly,  we  plotted  contours  of  the 
efficiencies  of  AV  and  SS  for  initial  numbers  of 
bombers  and  defenders.  From  these  plot  we  were  able  to 
see  how  the  initial  numbers  of  bombers  and  defenders  in 
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combat   affect   the   variability   of  uncertainty  in  the 

battle  outcomes  and  the  AV  and  SS  performance. 

a.   Variability  of  uncertainty  of  the  battle 
outcomes 

Equation     21     provides     an    adequate 

description  of  the  response   surface  that  characterizes 

the  variability  of  uncertain  in  the  expected  numbers  of 

live  defenders  at  the  end   of   the   battle   in   the  BCD 

model  .  The  response  ^Defender  ^s  the  expected  amount  of 

uncertainty  in  the  defender  estimate. 

vDefender=     -14-12    +    .<Z)J>6x<\     +    .  1  04x2    -     .009x5    -     .014x1**2 

+  .049x-,x2  +  .008x2**2  -  .051x1X3 
+  .040x2x5  -  .034x3**2  (21) 


TABLE  5.4  ANOVA  FOR  THE  EXPECTED  AMOUNT  OF  UNCERTAINTY 
IN  THE  DEFENDER  ESTIMATE 


SOURCE 


d.f 


SS 


MS 


F-RATIO 


Fitted  Surface 

9 

.  1302 

.0145 

8.06 

Lack  of  Fit 

3 

.01  02 

.0034 

1  .90 

Pure  Error 

2 

.0036 

.0018 

Total 

14 

.  1440 

R-Square= . 9043    Mean  Var iance= . 1 1 94    Std.  Dev.=.0525 


We  see  from  ANOVA  Table  5.4  that  the 
variation  in  the  variance  of  live  defenders  is 
insignificant  at  the  95%  level  (its  F-Ratio  of  8.06  is 
less  than  F    g^  .  g    2  =   19.38).  The   Lack  of   fit  is  also 
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insignificant  (F-Ratio=  1.90  <  F-g^.-j  2  =  19.16).  The 
R-Square  value  is  .9043  which  indicates  that  about  90% 
of  the  variation  in  the  variance  of  the  expected  number 
of  live  defenders  is  accounted  for  by  Equation  21 . 
Further  analysis  reveals  that  the  response  surface  is 
shaped  like  a  rising  ridge.  The  plots  of  contours  at 
Figure  5.1  illustrate  the  nature  of  this  Response 
surface.  These  pictures  show  that  initial  numbers  of 
bombers  and  defenders  affect  the  variance  of  defenders 
at  Time  25.  At  Times  75  and  125  the  initial  numbers  of 
bombers  have  little  influence.  Here  the  variance  of  the 
expected  number  of  live  defenders  is  affected  solely  by 
the  increase  in  the  number  of  initial  defenders.  Hence, 
as  the  number  of  defenders  increases,  the  variability 
of  uncertainty  in  the  estimate  of  the  expected  number 
of  live  defenders  at  the  end  of  the  battle  increases. 

Equation  22  provides  an  adequate 
description  of  the  response  surface  that  characterizes 
the  variability  of  uncertainty  in  the  expected  numbers 
of  live  bombers  at  the  end  of  the  battle  in  the  BCD 
model.  VBomber  is  the  expected  amount  of  uncertainty  in 
the  Bomber  estimate. 

vBomber=  •  1 483  +  -035x-|  -  .002x2  +  . 1 04x3  -  .031x1**2 

-  .041x1x2  +  .032x2**2  -  .067x-|X3 
+  .038x2x3  -  .007x3**2  (22) 
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Figure  5.1  Contour  Plots  of  the  Response  Surface  for 
the  Expected  Amount  of  Uncertainty  in  the  Defender 
Estimate 
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TABLE  5.5  ANOVA  FOR  THE  EXPECTED   AMOUNT  OF  UNCERTAINTY 
IN  THE  BOMBER  ESTIMATE 


SOURCE 


d.f . 


SS 


MS 


F-RATIO 


Fitted  Surface 
Lack  of  Fit 
Pure  Error 
Total 


9 

.  1331 

.0149 

6  .77 

3 

.0073 

.0024 

1  .08 

2 

.0045 

.0022 

4 

.  1408 

R-Square= . 9188    Mean  Var iance= . 1 1 88    Std.  Dev.=.0485 

ANOVA  Table  5.5  shows  that  the  variation  in 
the  variance  of  live  bombers  is  insignificant  at  the 
95%  level  (its  F-Ratio  of  6.77  is  less  than  F#g5.g  2  = 
19.38).  The  Lack  of  fit  is  also  insignificant  (F-Ratio= 
1.08  <  F. 95 «3,2  =  19.16).  The  R-Square  value  indicates 
that  about  92%  of  the  variation  in  the  variance  of  live 
bombers  is  accounted  for  by  Equation  22.  Further 
analysis  reveals  that  this  response  surface  is  also 
shaped  like  a  rising  ridge.  The  plots  of  contours  at 
Figure  5.2  illustrate  the  nature  of  this  response 
surface.  These  figures  show  that  the  variance  of  the 
expected  number  of  live  bombers  generally  increase  as 
the  initial  number  of  bombers  increases, 
b.   Antithetic  Variates. 

Equation  23  provides  an  adequate 
description  of  an  AV  response  surface  in  the  estimation 
of  the  expected  numbers  of  live  defenders. 
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VARIANCE  OF  BOMBER  ESTIMATE  AT  TIME  25 
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Figure  5.2  Contour  Plots  of  the  Response  Surface  for 
the  Expected  Amount  of  Uncertainty  in  the  Bomber 
Estimate 
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ADefender=  10-20  -  .09x-j  +  1.70x2-  .46x3  -  4.25x1**2 

+  .55x-|X2  -  1.68x2**2  -  .28x^3 
+  3.40x2x5  -  2.85x3**2  (23) 

We  see  from  ANOVA  Table  5.6  that  the  variation  in  the 
efficiency  of  AV  is  insignificant  at  the  95%  level  (its 
F-Ratio  of  1.46  is  less  than  F#95;g>2  =  19.38).  The 
Lack  of  fit  is  also  insignificant  (F-Ratio=.46  < 
F.95;3,2  =  19.16).  The  R-Square  value  is  .8497  which 
indicates  that  about  85%  of  the  variation  in  the 
Defender  RE  values  is  accounted  for  by  Equation  23. 
Here,  defender  is  the  expected  simulation  efficiency 
of  AV  generated  to  reduce  the  uncertainty  in  the 
Defender  estimate. 


TABLE  5.6  ANOVA  FOR  THE  EXPECTED  EFFICIENCY  OF  AV  (RE 
VALUE)  GENERATED  TO  REDUCE  THE  UNCERTAINTY  IN  THE 
DEFENDER  ESTIMATE 

SOURCE  d.f.        SS         MS        F-RATIO 

1  .46 
.47 


Fitted  Surface 

9 

168.31 

18.70 

Lack  of  Fit 

3 

4.08 

1  .36 

Pure  Error 

2 

25.69 

1  2.84 

Total 

14 

198.08 

R-Square= .8497    Mean  RE  Value=5.54    Std.  Dev.=2.44 

Further  analysis  reveals  that  the  response 
surface  is  shaped  like  a  hill  with  a  gentle  slope  on 
one  side  and  a  fairly  steep  slope  on  the  other  side. 
The   maximum   value   of   this  surface  occurs  in  the  BCD 
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scenario  that  begins  combat  with  50  defenders  and  40 
bombers  and  ends  the  battle  at  time  75.  The  plots  of 
contours  at  Figure  5.3  illustrate  the  nature  of  this  AV 
response  surface.  These  pictures  clearly  show  how  AV 
performs  in  different  scenarios  for  times  of  25,  75, 
and  125.  Beside  having  its  best  performance  in  a 
scenario  that  ends  at  time  75,  AV  appears  to  be  strong 
in  scenarios  that  initiate  the  air  battle  with  at  least 
30  defenders  and  30  bombers.  Its  weakest  performance 
seems  to  occur  in  those  scenarios  that  commence  combat 
with  no  more  than  30  defenders  and  30  bombers.  The 
plots  of  contours  show  that  the  efficiency  of  AV  is 
subversive  in  those  scenarios  whose  simulation 
initializes  the  air  battle  with  30  or  less  defenders 
and  40  or  more  bombers.  For  these  scenarios,  simulation 
efficiency  of  AV  may  often  increase,  instead  of 
decrease,  the  uncertainty  in  the  Defender  estimate.  We 
will  discuss  why  this  is  so  in  the  Experimental  Result 
Section.  Similar  analysis  of  the  AV  performance  is  made 
for  the  Bomber  RE  values.  An  adequate  description  of 
the  AV  response  surface  in  the  estimation  of  the 
bombers  is  characterized  by 

YBomber=  7.73  +  .44x-,  -  .69x2  +  2.05x3  -  2. 45x-|**2 
+  .18x-|X2  -  2.05x2**2  +  2.20x-|X3 
+  2.50x2x3  -  .53x3**2  (24) 
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Figure  5.3  Contours  Plots  of  the  Responses  Surface  for 
the  Expected  Efficiency  of  AV  Generated  to  Reduce  the 
Uncertainty  in  the  Defender  Estimate 
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YBomber  is  "the  expected  simulation  efficiency  of  AV  to 
reduce  the  uncertainty  in  the  Bomber  estimate.  The 
ANOVA  Table  5.7  indicates  that  the  proportion  of  the 
total  variation  in  the  Bomber  RE  values  accounted  for 
in  Equation  24  is  over  87%.  Furthermore,  this  variation 
is  insignificant  at  the  95%  level  (F-ratio  value  of 
8.19);  lack  of  fit  is  also  insignificant  (F-Ratio= 
1.61). 


TABLE  5.7  ANOVA  FOR  THE  EXPECTED  EFFICIENCY  OF  AV  (RE 
VALUE)  GENERATED  TO  REDUCE  THE  UNCERTAINTY  IN  THE 
BOMBER  ESTIMATE 


SOURCE 


d.f . 


SS 


MS 


F-RATIO 


Fitted  Surface 

9 

1 18.74 

13.19 

8.19 

Lack  of  Fit 

3 

14.17 

4.72 

1  .61 

Pure  Error 

2 

3.23 

1  .61 

Total 

14 

136. 14 

R-Square= .8722   Mean  RE  Value=5.05    Std.  Dev . = 1 . 87 

Examining  this  response  surface  further  we 
find  that  the  shape  of  the  surface  changes  over  time. 
The  plots  of  contours  depicted  in  Figure  5.4  show  that 
the  shape  of  the  surface  looks  like  a  saddle  at  Time 
25,  a  gentle  slope  at  Time  75,  and  a  uniformly  rising 
ridge  at  Time  125.  The  critical  values  for  this  surface 
also  change  as  its  shape  changes.  Most  notable  are  the 
values  for  maximum  efficiency.  At  Time  25,  maximum 
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EFFICIENCY  OF  AV  IN  BOMBER  ESTIMATE  AT  TIME  25 
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Figure  5.4  Contour  Plots  of  the  Response  Surface  for 
the  Expected  Efficiency  of  AV  Generated  to  Reduce  the 
Uncertainty  in  the  Bomber  Estimate 
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efficiency  occurs  in  scenarios  that  begin  fighting  with 
less  than  20  defenders  and  bomber.  By  Time  125,  maximum 
efficiency  has  shifted  to  the  scenarios  that  start 
with  at  least  40  defenders  and  bombers.  The  least 
amount  of  AV  reduction  occurs,  at  any  Time,  in  those 
scenarios  that  initialize  the  combat  simulation  with 
more  than  30  defenders  and  less  than  10  bombers. 

Here  is  a  summary  of  what  is  revealed  by 
the  above  analysis.  AV ,  in  general,  seems  to  be  the 
strongest  and  most  consistent,  and  equally-distributed 
between  closely-matched  pairs  of  bombers  and  defenders. 
Furthermore,  the  larger  the  evenly-matched  contest  the 
greater  the  variance  reduction.  When  the  defenders  and 
bombers  are  not  evenly  matched,  AV  is  not  as  consistent 
and  does  not  provide  equal  variance  reduction  in  the 
estimation  of  the  pair  of  parameters.  It  is  strong  in 
the  estimation  of  the  larger  combatants  and  weak  in  the 
estimation  of  the  smaller  ones, 
c.   Stratified  Sampling. 

We  analyze  the  efficiency  of  SS  in  the 
simulation  of  the  BCD  model  in  a  similar  manner  as  we 
analyzed  the  efficiency  of  AV .  If  we  analyze 
Equations  25  and  26  in  terms  of  ANOVA  Tables  5.8  and 
5.9  respectively,  we  will  get  results  similar  to  those 
we  obtained  in  the  last  section.  Therefore,  we  will 
forego  this  particular  analysis.  Note  that  ^j)efender    is 
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the  expected  simulation  efficiency  of  SS  to  reduce  the 
uncertainty  in  the  Defender  estimate,  and  YBomber  is 
the  expected  simulation  efficiency  of  SS  to  decrease 
the  uncertainty  in  the  Bomber  estimate. 


YDefender=  5. 03  -  -39X1  +  .49x2  +  .08x3  -  .15x1»*2 
+  .22x-|X2  -  .70x2**2  -  . 45x^3  +  . 55x2x5 
-  .58x3**2  (25) 


TABLE  5.8  ANOVA  FOR  THE  EXPECTED  EFFICIENCY  OF  SS  (RE 
VALUE)  GENERATED  TO  REDUCE  THE  UNCERTAINTY  IN  THE 
DEFENDER  ESTIMATE 


SOURCE 


d.f 


SS 


MS 


F-RATIO 


Fitted  Surface 

9 

8.24 

.92 

2.49 

Lack  of  Fit 

3 

.53 

.  18 

.47 

Pure  Error 

2 

.75 

.37 

Total 

14 

9.52 

R-Square= .86661 


Mean  RE  Value=2.27 


Std.  Dev.=.51 


YBomber=  2  •  90  -  .06^  -  -03x2  +  .21x3  -  .41x^*2 

-  .20x-|X2  -  .64x2**2  +  .58x-,X3  +  .  50x2x5 

-  .36x3**2  (26) 
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TABLE  5-9  ANOVA  FOR  THE  EXPECTED  EFFICIENCY  OF  SS  (RE 
VALUE)  GENERATED  TO  REDUCE  THE  UNCERTAINTY  IN  THE 
BOMBER  ESTIMATE 


SOURCE 


d.f  . 


SS 


MS 


F-RATIO 


Fitted  Surface  9 

Lack  of  Fit  3 

Pure  Error  2 

Total  14 


5.18 

.58 

19.33 

1  .82 

.61 

20.33 

.06 

.03 

7.06 

R-Square= .7340    Mean  RE  Value=2.15     Std.  Dev.=.61 

The  contour  plots  at  Figures  5.5  and  5.6  appear  to 
have  similar  features.  They  show  a  relatively  flat 
surface  except  at  the  corners.  The  corner  with  50 
Bombers  and  50  Defenders  has  the  highest  response  and 
the  other  corners  have  low  response.  These  plots 
suggest  that  the  SS  performance  is  generally  consistent 
in  all  but  a  few  scenarios  in  the  BCD  model.  Maximum 
efficiency  of  SS  occurs  in  those  scenarios  that 
initialize  simulation  with  equally  large  numbers  of 
bombers  and  defenders.  It  is  very  weak  in  those 
scenarios  that  begin  combat  with  either  less  than  10 
defenders  and  more  than  40  bombers  or  more  than  40 
defenders  and  less  than  10  bombers. 
3 •   Experimental  Results 

The  experimental  results  can  be  summarized  in 
Tables  5.10  and  5.11.  Table  5.10  shows  the 
relationships  between  the  AV  performance  and  the 
uncertainty  in  the  Defender  estimate  and  the  SS 
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EFFICIENCY  OF  SS  IN  DEFENDER  ESTIMATE  AT  TIME  25 
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Figure  5.5  Contour  Plots  of  the  Response  Surface  for 
the  Expected  Efficiency  of  SS  Generated  to  Reduce  the 
Uncertainty  in  the  Defender  estimate 
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EFFICIENCY  OF  SS  IN  BOMBER  ESTIMATE  AT  TIME  25 
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Figure  5.6  Contour  Plots  of  the  Response  Surface  for 
the  Expected  Efficiency  of  SS  Generated  to  Reduce  the 
Uncertainty  in  the  Bomber  Estimate 
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performance  and  the  uncertainty  in  the  Defender 
estimate.  Similarly,  Table  5.11  shows  the  relationships 
between  the  AV  performance  and  the  uncertainty  in  the 
Bomber  estimate  and  the  SS  performance  and  the 
uncertainty  in  the  Bomber  estimate. 

We  further  examine  this  relationship  by 
analyzing  the  data  that  measure  the  uncertainty  (crude 
variance)  and  appropriate  variance  reduction  (RE 
Values)  in  Appendix  E.  After  we  applied  a  logarithmic 
transformation  to  the  data,  we  regress  the  RE  values 
on  the  crude  variance  data  and  observe  a  strong 
logarithmic  linear  relationship  between  uncertainty  and 
variance  reduction. 


TABLE  5.10  RELATIONSHIP  BETWEEN  UNCERTAINTY  AND  THE 
EFFICIENCY   OF   VARIANCE   REDUCTION    IN   THE   DEFENDER 
ESTIMATE 


INITIAL 

UNCERTAINTY 

VARIANCE 

REDUCTION 

Defenders 

Bombers 

Variance 

AV 

SS 

10 

10 

medium 

strong 

fair 

10 

30 

medium 

fair 

fair 

10 

50 

small 

weak 

weak 

30 

10 

large 

strong 

fair 

30 

30 

large 

strong 

strong 

30 

50 

medium 

fair 

fair 

50 

10 

large 

strong 

weak 

50 

30 

large 

strong 

fair 

50 

50 

large 

strong 

fair 
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TABLE  5.11   RELATIONSHIP  BETWEEN  UNCERTAINTY  AND  THE 
EFFICIENCY  OF  VARIANCE  REDUCTION  IN  THE  BOMBER  ESTIMATE 


INITIAL 

UNCERTAINTY 

VARIANCE 

REDUCTION 

Defenders  B 

ombers 

Variance 

AV 

SS 

10 

10 

medium 

strong 

fair 

10 

30 

large 

strong 

fair 

10 

50 

large 

strong 

fair 

30 

10 

smal  1 

fair 

fair 

30 

30 

large 

strong 

fair 

30 

50 

large 

strong 

fair 

50 

10 

smal  1 

weak 

weak 

50 

30 

medium 

fair 

fair 

50 

50 

large 

strong 

fair 

This  relationship  is  manifested  in  the 
multiplicative  equation  shown  in  Table  5.12.  VAR  is  the 
value  of  uncertainty  or  variance  of  the  corresponding 
estimate  obtained  from  crude  simulation,  and  Y  is  the 
simulation  efficiency  of  the  variance  reduction  or  RE 
value  for  the  corresponding  estimate. 


TABLE  5.12  ANALYSIS  OF  THE  RELATIONSHIP  BETWEEN 
UNCERTAINTY  AND  EFFICIENCY  OF  VARIANCE  REDUCTION  IN 
PARAMETER  ESTIMATES 

CORRELATION     SET 
VRT   ESTIMATE      RELATIONSHIP       COEFFICIENT    POINT 

AV  Defender  Y=  10.43  *  VAR**. 362 

AV  Bomber  Y=  9.73  *  VAR**. 362 

SS  Defender  Y=  3.18  *  VAR**. 144 

SS  Bomber  Y=  2.99  ■  *  VAR**. 147 


.8609 

.00154 

.8217 

.00186 

.5601 

. 00032 

.6054 

. 00058 
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The  correlation  coefficient  reveals  the 
strength  of  the  logarithmic  linear  relationships 
between  uncertainty  and  the  efficiencies  of  AV  and  SS . 
With  the  values  of  the  exponent  in  the  equations  being 
less  than  one  and  the  values  of  the  correlation 
coefficient  being  positive,  the  efficiencies  of  AV  and 
SS  are  observed  to  increase,  at  a  decreasing  rate,  as 
the  uncertainty  (variance)  increases.  We  obtained  the 
Set  Point  by  setting  Y=1  in  the  corresponding  equation 
and  solving  for  VAR .  At  this  value,  simulation 
efficiency  nether  increases  nor  decreases.  Now  if  we 
observe  a  value  of  uncertainty,  or  variance  obtained 
from  crude  simulation,  above  this  set  point,  then  we 
expected  to  get  an  efficiency  of  a  VRT  to  increase  the 
simulation  efficiency.  On  the  other  hand,  if  the  value 
is  below  the  set  point,  then  we  expect  the  efficiency 
of  the  VRT  to  decrease  the  simulation  efficiency. 

Here   is    the   bottom    line   on    AV   and   SS 

performance  in  the  BCD  model; 

1 .   If  we  apply  antithetic  variates  to  the  simulation 
of  the  BCD  model , 

a.   We  may  expect  the  variability  of   uncertainty 
in  the  defender  estimate 

( 1 )   to  decrease  if  the  variance  of  the 

estimate  obtained   from  crude  simulation 
is  at  least  .00154,  and 
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(2)   to  decrease,  at  a  decreasing   rate,  with 
an  increase  in  the  variance  of  the 
estimate  obtained   from  the  crude 
simulation . 

b.   We  may  expect  the  variability  of   uncertainty 
in  the  bomber  estimate 

( 1 )  to  decrease  if   the  variance   of  the 
estimate  obtained  from  the  crude 
simulation  is  at  least  .00186,  and 

(2)  to  decrease,  at  a  decreasing  rate,  with 
an  increase  in  the  variance  of  the 
estimate  obtained  from  the  crude 
simulation . 

2.   If  we  apply  stratified  sampling  to  the  simulation 
of  the  BCD  model , 

a.  We  may  expect  the  variability  of   uncertainty 
in  the  defender  estimate 

(1  )   to  decrease  if  the  variance  of  the 

estimate  obtained   from  crude  simulation 
is  at  least  .00032,  and 

(2)   to  decrease,  at  a   decreasing  rate,  with 
an  increase   in  the  variance  of  the 
estimate  obtained  from  the  crude 
simulation . 

b.  We  may  expect  the  variability  of   uncertainty 
in  the  bomber  estimate 

(1 )  to  decrease  if  the  variance  of  the 
estimate  obtained  from  the  crude 
simulation  is  at  least  .00058,  and 

(2)  to  decrease,  at  a   decreasing  rate,  with 
an  increase  in  the  variance  of  the 
estimate  obtained  from  the  crude 
simulation . 

E .   SUMMARY 

We   illustrated   the   performance  of  the  antithetic 

variates  and  stratified  sampling   in  the   simulation  of 

the   BCD    model.   The    manifestation   of   their   pair 
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performance  was  characterized  by  response  surface 
equations  and  plots  of  contour  lines.  Both  VRTs  were 
shown  to  be  effective  in  increasing  simulation 
efficiency,  but  they  perform  differently  in  the  BCD 
model.  AV  provides  the  largest  amount  of  variance 
reduction  but  is  more  volatile.  AV  increases  the 
simulation  efficiency  on  the  average  of  5  times  the 
crude  simulation;  it  is  strong  in  the  BCD  scenarios 
where  there  is  large  amount  of  uncertainty  in  the 
battle  outcomes  for  live  bombers  and  defenders,  and 
weak  in  those  scenarios  where  there  is  little  amount  of 
uncertainty  in  the  battle  outcomes.  SS ,  on  the  hand, 
has  a  more  consistent  performance.  SS  increases  the 
simulation  efficiency  at  a  mean  of  2  times  the  crude 
simulation.  It  performs  nearly  the  same  in  every 
scenario  except  where  the  uncertainty  is  large. 
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VI .   CONCLUSIONS  AND  RECOMMENDATIONS 

A.   CONCLUSIONS 

The  objective  of  this  pilot  study  has  been  to 
investigate  the  applicability  of  VRTs  to  reduce  the 
inherent  variability  in  stochastic  combat  models.  We 
examined  the  effects  of  applying  AV  and  SS  to  the 
simulation  of  a  simple  stochastic  combat  model.  We  have 
now  shown  that  AV  and  SS  are  applicable  to  this 
stochastic  combat  model.  We  can  infer  that  these  VRTs 
are  indeed  capable  of  working  in  stochastic  combat 
models,  and  their  prospects  in  larger  and  more  complex 
stochastic  combat  models  are  even  more  promising.  The 
conditions  of  monotonicity  and  synchronization  are 
essential  parts  of  the  design  of  the  simulation  program 
for  these  models.  Hence,  we  feel  that  sizable  increase 
in  simulation  efficiency  is  possible  if  these 
requirements  are  met  in  the  simulation. 

The  experimental  results  of  applying  VRTs  to  the 
BCD  simulation  model  show  that  the  strength  of  AV  and 
SS  is  influenced  by  uncertainty.  A  strong  variance 
reduction  results  from  a  large  variance  of  the  estimate 
obtained  from  crude  simulation.  A  weak  variance 
reduction  is  caused  from  a  small  variance  of  the 
estimate  obtained  from  crude  simulation.  Hence,  sizable 
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and  consistent  variance  reduction  depends  on  large 
variability  of  the  simulated  output  from  the  stochastic 
combat  model.  Therefore,  the  variability  of  the  output 
from  larger  and  more  complex  stochastic  models  must 
also  be  large  enough  to  obtain  the  size  and  consistency 
of  simulation  efficiency  and  variance  reduction  one 
desires  from  the  applications  of  these  VRTs  to  such 
models . 

B .   RECOMMENDAT IONS 

The  pilot  study  presented  in  this  thesis  provides 
a  base  for  further  studies  in  the  applications  of  AV 
and  SS  to  large-scale,  real  world,  stochastic  combat 
simulation  models.  Usually  complex  simulation  models 
have  many  subroutines  or  modules.  The  variability  of 
uncertainty  in  the  output  data  from  these  modules  may 
vary  from  low  to  high.  We  recommend  that  a  study  of 
this  matter  focus  on  the  degree  of  variability  of 
uncertainty  in  the  output  data  from  each  module.  The 
interest  of  the  study  should  be  concerned  with  the 
relationship  between  the  performance  of  the  VRT  and  the 
variability  of  the  output  data  from  each  module.  The 
results  should  indicate  where  and  how  the  VRT  may  be 
used  in  the  model  in  order  to  maximize  the  simulation 
efficiency  of  the  model.  For  example,  if  the  study 
shows   that   a   VRT   performs   strongly  in  a  particular 
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module  whose  variability  of  output  data  is  large  and  it 
performs  poorly  in  another  module  whose  variability  of 
output  is  small,  then  the  study  should  recommend  that 
the  VRT  be  used  in  the  module  in  which  it  performs  best 
and  not  be  used  in  that  module  for  which  it  performs 
poorly.  Using  VRTs  in  the  module  with  small  variability 
would  most  likely  decrease  simulation  efficiency  for 
that  model  and,  at  worst,  suboptimize  the  overall 
performance  of  the  VRT  in  the  complex  combat  model . 
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APPENDIX  A- 
FORMULATION  OF  THE  BCD  MARKOVIAN  MODEL 

by 

Professor  Donald  P.  Gaver 
Naval  Postgraduate  School 
Monterey,  California 

B  bombers  are  approaching  a  group  of  D  defenders. 
When  the  two  groups  approach  within  range  each  defender 
searches  for  a  bomber;  after  he  finds  one  they  engage 
In  combat.  Either  bomber  or  defender  may  win  the 
combat;  the  survivor  becomes  "free",  and  is  a  candidate 
for  the  next  combat.  In  general,  bombers  attempt  to 
avoid  combat,  defenders  seek  it  out. 

This  situation  becomes  a  tri-variate  Markov  chain 
if  the  following  state  is  defined:  { B( t  )  , C( t ) , D( t ) } . 
Here  t  is  conveniently  measured  from  the  time  bombers 
and  defenders  are  close  enough  to  permit  combat  at  all, 
B(t)  is  the  number  of  free  bombers  at  time  t 
thereafter;  ditto  for  D(t),  the  number  of  free 
defenders;  C(t)  is  the  number  of  one-on-one  combats. 
Here  are  a  set  of  transition  rates: 

(1)  Combats  begin.  If  B(t)=b,  C(t)=c,  D(t)=d  and  if  f 
is  the  rate  at  which  free  bombers  are  found  by 
free  defenders  then  with  probability 
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((e-jt)b)d  =  1  -  fbdt  -+  o(t) 

no  defender   finds  a  bomber  in  time  t.  Hence  the 
probability  that  a  defender  does  find  a  bomber  is 
ibdt  +   o(t).  This  is  the  rate  at  which  free 
bombers  and  defenders  get  converted  to   combats: 
the  state 

(b,  c,  d)  ->  (b-1  ,  c  +  1  ,  d-1)   with  prob  ibdt. 

(2)  Defenders   win.  Same  initial  conditions.  If  a 
bomber  is  in  combat  with  a  defender  the 
probability  that   a  defender  shoots  down  the 
bomber  in  time  t  (combat  duration)   if  the  latter 
doesn't  hit  the  defender  is 

P(Combat  duration  i    t  !  bomber  doesn't  hit) 
=  1  -  e_at. 

Likewise,  the   probability  that  the  bomber  shoots 
down  the   defender  is 

P{ Combat  duration  i    t  !  defender  doesn't  shoot) 
=  1  -  e-0t. 
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Now  suppose  they  both  shoot »  doing  so 
independently.  Model   the  probability  that  both 
survive  to  t  as 

P<Combat  duration  >  t)  =  (  e~at )  (  e"Bt )  . 

Now  the  probability  that  combat  lasts  until  t  and 
is  terminated   by  a  defender  shooting  down  the 
bomber  is 

P{Combat  ends  is  (  dt  )  ,  Defender  wins} 

=  (e_atadt)e"et 

=  (e-(«  +  n)t(0  +  ,))(„/(„+,). 

This  shows  that  a  single  combat  duration  is 
exp(a  +  9)  and  the  event  of  a  defender's  winning 
is  independently  a/(a  +  e).  Likewise,  the  combat 
duration  is   exp(a  +  e )  and  a  bomber's  win  is, 
independently,  e/(  +  e).  If  there  are  c  combats 
going  on   then  the  first  combat  to  ends  does  so 
in  time  exp(c(a  +  0)). 
Hence 

(b,  c,  d)   *  (b,  c-1 ,  d+1 )  with  prob  acdt 
(b,  c,  d)  *  (b+1 ,c-1 ,  d)  with  prob  cdt. 
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(3)   Simulation  and  Sojourn.  The  above  shows   that  we 

may  simulate  the  combat  as  follows. 

(i)      You  are  in  state  (b,  c,  d ) .  Obtain  a 
sojourn  time  in  that  state  that  is 

exp($bd  +  (a  +  e)c) 

i.e.    Sbcd  =  1  /  (ibd  +  (a  +  r>)c). 

The  system  stays  in  state  (b,c,d)   for 
time  [0,Sbcd). 

(ii)    With  probability 

ibd  /  (ibd  +  (a  +  b)c)  ♦  (b-1  ,  c  +  1  ,  d-1) 

(NEW  COMBAT)  at  time  Sbcd. 
(iii)   With  probability 

aC  /   (ibd  +  (a  +  e)c)  *    (b,   C-1,   d+1  ) 

(DEFENDER  SHOOTS  DOWN  BOMBER)  at 
time  Sbcd. 

(iv)    With  probability 

ec  /   (ibd  +  (a  +  p)c)  *  (b  +  1,c-1,d) 


(BOMBER  SHOOTS  DOWN  DEFENDER)  at 
time  Sbcd. 
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APPENDIX  B 

FORTRAN  PROGRAM  LISTING  FOR  THE  CRUDE  SIMULATION 

OF  THE  BCD   MODEL 


DIMENSION  BX( 101 ) ,DX( 101 ) ,SEED(2) ,B0X(2) ,D0X(2) 

INTEGER  I,N,BB,DD,R 

REAL*4  X,Y,Z,BX,DX,TXT,BOX,DOX 

DOUBLE  PRECISION  SEED 

DATA  SEED  / 1 234 . 0D0 , 567890 1 23 . 0D0/ 

R=  2 

N=1  00 
C 

C   RECEIVE  INPUT  DATA  FROM  TERMINAL 
C 

WRITE( * ,3) 

3  FORMAT ( 1 X, 'ENTER  THE  NUMBER  OF  BOMBERS') 
READ  '  ( 12)  '  ,  BB 

WRITE( * ,4) 

4  FORMAT ( 1X, 'ENTER  THE  RATE  WHICH  A   BOMBER  SHOOTS 
DOWN  A  DEFENDER ' ) 

READ  ' (F5.3) ' ,  Y 
WRITE( *  ,5) 

5  FORMAT ( 1X, 'ENTER  THE  NUMBER  OF  DEFENDERS') 
READ  '  (  12)  '  ,  DD 

WRITE( * ,6) 

6  FORMAT ( 1X, 'ENTER  THE  RATE  WHICH  A  DEFENDER' 
&  'SHOOTS  DOWN  A  BOMBER') 

READ  ' (F5.3)  '  ,  X 

WRITE( * ,7) 'ENTER  THE  RATE  WHICH  FREE  DEFENDERS' 
&  'FIND  FREE  BOMBERS' 

7  FORMAT (1X, A) 
READ  ' ( F5 . 3  )  '  ,  Z 

WRITE( * ,8) 'ENTER  THE   TIME  DURATION  OF  THE 
&  'BATTLE' 

8  FORMAT (1X, A) 

READ  ' (F5.3)  '  ,  TXT 
C 

C    RUN  REPLICATIONS  OF  N  BATTLES  AND  OBTAIN  SUMMARY  OF 
C    N  BATTLES 
C 

CALL  BATTLE( N , R , SEED , X , Y , Z , BB , DD , TXT , BX , DX ) 
C 

C   COMPUTE  STATISTICAL  OUTPUT  ANALYSIS  OBTAIN  PARAMETER 
c   ESTIMATES 
C 

CALL  STAT(N,R,BX,DX,BOX,DOX) 


86 


C   FORMAT  AND  PRINT  OUTPUT  OF  PARAMETER  ESTIMATES 
C 

WRITE(3,279)  N 

279  FORMAT (15X,  'SAMPLE  SIZE '  ,  16  ,  / 15X  ,  6 (  ' - '  )  , IX  , 
&   4(  '-'  )) 

WRITE(3,280)  BB,DD 

280  FORMAT( //1X.I4 ,2X, 'BOMBERS' , 6X , 'VERSUS' , 6X , 
&  14 ,2X,  'DEFENDERS'  ,/lX,13(  '-'  )  , 18X,14(  '-'  )  ) 

290    WRITE  (3,300)TXT 

30  0    FORMAT ( //18X, 'TIME' ,F6.1,/,47('-')) 

WRITE(3,310 ) 
310    FORMAT(19X, 'BOMBER' , 2X , 'DEFENDER' ,/18X, 
&  7(  '  =  '  ),2X,8(  '  =  '  )) 

WRITE(3,320)  BOX ( 1  )  , DOX ( 1 ) 
320    FORMAT( IX, 'AVERAGE' ,7X,2F10.4 ) 

WRITE(3,330)  BOX( 2  )  , DOX( 2 ) 
330    FORMAT ( IX, 'VARIANCE' ,6X,2F10.4 ) 
STOP 
END 
C 
C 

C    SUBROUTINE  BATTLE 
C 

SUBROUTINE 
&  BATTLE ( N , R , SEED , X , Y , Z , BB , DD , TXT , BX , DX ) 
INTEGER  BB,DD,I ,K,N,R 

REAL  X ( N+l ) , DX ( N+l ) , GA , SO J ,X,Y,Z,B,C,D,NC, 
&  BW,DW, INF, T, TIME, TXT 

DOUBLE  PRECISION  SEED(R) 
C 

C   INITIALIZE  STATISTICAL  COUNTERS 
C 

BX(N+1)=  0.0 
DX(N+1)=  0.0 
C 

C   RUN  N  REPLICATIONS 
C 

DO  200  I-1,N 
C 

C    INITIALIZE  START-TO-BATTLE  VALUES 
C 

T-0.0 
C-0.0 

B=REAL(BB) 
D=REAL(DD) 
BW=  0.0 
DW=  0.0 
NC=  1.0 
C 

C   OBSERVE  OCCURRENCE  OF  AN  EVENT 
C 
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100  CALL  UNIFOR(SEED( \ ) ,  GA , 1 ) 

C 

C  DETERMINE   NEXT  INTERIM  EVENT  AND  UPDATE  THE   STATE 
C   (B,C,D) 
C 

IF(GA  .LE.  BW)  THEN 
B=B  +1.0 
C=C  -1.0 
D  =  D 
ELSE  IF  (GA  .LE.  ( BW+NC ) )  THEN 
B=B  -1.0 
C  =  C  +1.0 
D=D  -1.0 
ELSE 
B  =  B 

C=C  -1.0 
D=D  +1.0 
END  IF 
C 

C   COMPUTE  MEAN  TIME  IN  STATE  (B,C,D) 
C 

IF  ((B  .EQ.  0.0  .OR.  D  .EQ.  0.0)   .AND.  C 
&         .EQ.  0.0)  THEN 

INF=  1000000.0 
ELSE 

INF=  1.0/(Z*D*B  +  (X  +  Y)*C) 
END  IF 
C 

C   GENERATE  SOJOURN  TIME  IN  STATE  (B,C,D) 
C 

CALL  EXP0N(SEED(2) , SO J , 1 ) 
TIME=  -INF  *  ALOG(SOJ) 
C 

C   ADVANCE  THE  SIMULATED  TIME  OF  THE  AIR  BATTLE 
C 

T=  T  +  TIME 
C 
C   COMPUTE  PROBABILITY  OF  NEXT  INTERIM  EVENTS  OCCURRING 


BW=  Y*C*INF 

DW=  X*C*INF 

NC=  Z*B*D*INF 
C 
C   CHECK  CONDITIONS  FOR  OCCURRENCE  OF  TERMINATION  EVENT 


IF  (T  .LT.  TXT)  GOTO  100 
C 

C    ACCUMULATE  SUMMARY  OF  N  BATTLE  OUTCOMES 
C 
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BX( I)=  B  +  C 
DX( I)=  D  +  C 
BX(N+1 )=  BX(N+1 )  +  BX(I) 
DX(N+1 )=  DX(N+1 )  +  DX( I) 
200    CONTINUE 

RETURN 

END 
C 
C 

C    SUBROUTINE  EXPON 
C 

SUBROUTINE  EXPON( SEED2 , A2 , K ) 

INTEGER  I,K 

REAL  STIM.A2, INF 

DOUBLE  PRECISION  EFF,SEED2 

EFF=  21 47483647. 0D0 

SEED2=DMOD( 1 6807.0D0  *  SEED2,EFF) 

A2=  SEED2/EFF 

RETURN 

END 
C 
C 

C    SUBROUTINE  UNIFOR 
C 

SUBROUTINE  UNIFOR( SEED1 , A1 , K ) 

INTEGER  I,K 

REAL  A1 

DOUBLE  PRECISION  EFF,SEED1 

EFF=  2147483647. 0D0 

SEED1 =DMOD( 1 6807.0D0  *  SEED1,EFF) 

A1=  SEED1/EFF 

RETURN 

END 

C 

C 

C    SUBROUTINE  STAT 

C 

SUBROUT INE  STAT ( N , R , BX , DX , BOX , DOX ) 
INTEGER  J,R,N 

REAL  BX(N+1  ),DX(N+1  )  , BOX ( R ) , DOX ( R ) 
B0X(2)=  0.0 
DOX(2)=  0.0 
C 

C    COMPUTE  THE  ESTIMATES  OF  THE  SAMPLE  MEAN  AND 
C   VARIANCE 
C 

BOX( 1 )  =  BX(N+1 )/N 
DOX( 1 )  =  DX(N+1 )/N 
DO  260  1=1 ,N 

B0X(2)=  BOX(2)  +  (BX( I)-BOX( 1 ) )**2 
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D0X(2)=  D0X(2)  +  (DX(I)-DOX( 1 ))**2 
260       CONTINUE 

B0X(2)=  B0X(2)/(N*(N-1 ) ) 
D0X(2)=  DOX(2)/(N*(N-1 ) ) 
RETURN 
END 
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APPENDIX  C" 

FORTRAN  PROGRAM  LISTING  OF  THE  BCD  SIMULATION 
USING  ANTITHETIC  VARIATES 


DIMENSION  BX( 101  )  ,DX( 101  ) ,SEED(2) ,B0X(2)  ,DOX(2) 

INTEGER  I,N,BB,DD,R 

REAL*4  X,Y,Z,BX,DX,TXT,BOX ,DOX 

DOUBLE  PRECISION  SEED 

DATA  SEED  / 1 234 . 0D0 , 5678901 23 . 0D0/ 

R=  2 

N  =  50 
C 

C   RECEIVE  INPUT  DATA  FROM  TERMINAL 
C 

WRITE( * ,3) 

3  FORMAT ( 1 X, 'ENTER  THE  NUMBER  OF  BOMBERS') 
READ  ' ( 12)  '  ,  BB 

WRITE( *  ,4) 

4  FORMAT ( 1X, 'ENTER  THE  RATE  WHICH  A  BOMBER  SHOOTS' 
&  'DOWN  A   DEFENDER') 

READ  ' (F5.3)  '  ,  Y 
WRITE( * ,5) 

5  FORMAT ( 1X, 'ENTER  THE  NUMBER  OF  DEFENDERS') 
READ  ' ( 12)  '  ,  DD 

WRITE( * ,6) 

6  FORMAT ( 1X, 'ENTER  THE  RATE  WHICH  A  DEFENDER 
&  SHOOTS  DOWN  A  BOMBER') 

READ  ' (F5.3)  '  ,  X 

WRITE( * ,7) 'ENTER  THE  RATE  WHICH   FREE  DEFENDERS' 
&  'FIND  FREE  BOMBERS' 

7  FORMAT (1X, A) 
READ  ' (F5.3) ' ,  Z 

WRITE( * ,8) 'ENTER  THE   TIME  DURATION  OF  THE' 
&  'BATTLE' 

8  FORMAT (1X, A) 

READ  ' (F5.3) ' ,  TXT 
C 

C    RUN  REPLICATIONS  OF  N  BATTLES  AND  OBTAIN  SUMMARY  OF 
C    N    BATTLES 
C 

CALL  BATTLE( N , R , SEED , X , Y , Z , BB , DD , TXT , BX , DX ) 
C 

C   COMPUTE  STATISTICAL  OUTPUT  ANALYSIS  OBTAIN  PARAMETER 
C   ESTIMATES 
C 

CALL  STAT(N,R,BX,DX,BOX,DOX) 
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C   FORMAT  AND  PRINT  OUTPUT  OF  PARAMETER  ESTIMATES 
C 

WRITE(3,279)  N 

279  FORMAT ( 1 5X, 'SAMPLE  SIZE '  , 16  ,  / 1 5X , 6(  ' - '  )  , 
&  1X,4(  '-'  )) 

WRITE(3,280)  BB,DD 

280  FORMAT(//1X, 14, 2X, 'BOMBERS' , 6X , 'VERSUS' , 

&  6X, I4,2X,  'DEFENDERS'  , /1X, 13(  '-'  ),18X,14(  '-'  )) 
290    WRITE  (3,300)TXT 
300    FORMAT ( //18X,  'TIME'  , F6 . 1  , / , 47 (  ' - '  )  ) 

WRITE(3,310) 
31 0    FORMAT ( 1 9X, 'BOMBER' ,  2X, 'DEFENDER' , 
&  /18X,7(  '  =  '  ) ,2X,8(  '  =  '  )  ) 

WRITE(3,320)  BOX( 1  )  ,DOX( 1  ) 
320    FORMAT( 1X, 'AVERAGE' ,7X,2F10.4) 

WRITE(3,330)  BOX(2)  ,D0X(2) 
330    FORMAT( 1 X, 'VARIANCE' ,6X,2F1 0.4) 
STOP 
END 
C 
C 

C    SUBROUTINE  BATTLE 
C 

SUBROUTINE 
&  BATTLE( N , R , SEED ,X,Y,Z,BB,DD, TXT , BX , DX ) 
INTEGER  BB,DD,H,I,J,W,K,N,R 
REAL  GA(2, 1 000) , S0J(2, 1 000) , BX(N+1 ) , 
&  DX(N+1 ) ,BAT(50,2) ,DAT(50,2) , 
&  X,Y,Z,T, TXT, TIME, INF , NC , BW , DW , B , C , D 
DOUBLE  PRECISION  SEED(R) 
K=  1000 
C 

C  INITIALIZE  STATISTICAL  COUNTERS 
C 

BX(N+1 )=  0.0 
DX(N+1  )=  0.0 
C 

C  RUN  N  REPLICATIONS 
C 

DO  200  1=1 ,N 

CALL  SOJOUR(SEED( 1 ) ,SOJ,R,K) 
CALL  STATE(SEED(2) ,GA,R,K) 
DO  175  J=1 ,R 
C 

C    INITIALIZE  START-TO-BATTLE  VALUES 
C 

H  =  0 

T  =  0  .  0 

C  =  0.0 

B=REAL(BB) 

D=REAL(DD) 
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BW=  0.0 
DW=  0.0 
NC=  1 .0 
C 

C   OBSERVE  OCCURRENCE  OF  NEXT  EVENT 
C 

100  H  =  H+1 

C 

C  DETERMINE  NEXT  INTERIM  EVENT  AND  UPDATE  THE  STATE 
C   (B,C,D) 
C 

IF(GA(J,H)  .LE.  BW)  THEN 
B=B  +1.0 
C=C  -1.0 
D  =  D 
ELSE  IF  (GA(J,H)  . LE .  ( BW+NC ) )  THEN 
B=B  -1.0 
C=C  +1.0 
D=D  -1.0 
ELSE 
B  =  B 

C=C  -1.0 
D=D  +1.0 
END  IF 
C 

C   COMPUTE  MEAN  TIME  IN  STATE  (B,C,D) 
C 

IF  ((B   .EQ.  0.0  .OR.  D  .EQ.  0.0)  .AND.  C 
&  .EQ.  0.0)  THEN 

INF=  1000000.0 
ELSE 

INF=  1.0/(Z*D*B  +  (X  +  Y)*C) 
END  IF 
C 

C   COMPUTE  SOJOURN  TIME  IN  STATE  (B,C,D) 
C 

TIME=  -INF  *  ALOG(SOJ( J,H)  ) 
C 

C   ADVANCE  THE  SIMULATED  TIME  OF  THE  AIR  BATTLE 
C 

T=  T  +  TIME 
C 
C   COMPUTE  PROBABILITY  OF  NEXT  INTERIM  EVENTS  OCCURRING 


BW=  Y*C*INF 

DW=  X*C*INF 

NC=  Z*B*D*INF 
C 

C   CHECK  FOR  OCCURRENCE  OF  TERMINATION  EVENT 
C 
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IF  (T  .LT.  TXT)  GOTO  100 
C 

C    RECORD  RESULTS  OF  BATTLE 
C 

BAT( I, J)=  B  +  C 
DAT( I , J)=  D  +  C 
175       CONTINUE 
C 

C    ACCUMULATE  SUMMARY  OF  N  BATTLE  OUTCOMES 
C 

BX( I  )  =  (BAT( I , 1  )  +  BAT (1,2))*. 5 
DX( I)=(DAT( I, 1 )  +  DAT (1,2))*. 5 
BX(N+1 )=  BX(N+1)  +  BX(I,J) 
DX(N+1 )=  DX(N+1 )  +  DX(I,J) 
200    CONTINUE 
RETURN 
END 
C 
C 

C    SUBROUTINE  SOJOUR 
C 

SUBROUTINE  SOJOUR( SEED2 , A2 , W , K ) 

INTEGER  I,W,K 

REAL  A2(W,K) 

DOUBLE  PRECISION  EFF , SEED2 

EFF=  21 47483647. 0D0 

DO  10  1  =  1  ,K 

SEED2=DMOD( 1 6807.0D0  *  SEED2,EFF) 
A2( 1 , I )=  SEED2/EFF 
A2(2,  I)=  1.0  -  A2( 1  , I) 
10     CONTINUE 
RETURN 
END 
C 
C 

C    SUBROUTINE  STATE 
C 

SUBROUTINE  STATE( SEED1 , A1 , W , K ) 

INTEGER  I,K,W 

REAL  A1 (W,K) 

DOUBLE  PRECISION  EFF.SEED1 

EFF=  21 47483647. 0D0 

DO  10  1  =  1  ,K 

SEED1 =DMOD( 1 6807.0D0  *  SEED1.EFF) 
A1 ( 1 ,I)=  SEED1 /EFF 
A1 (2,1)=  1.0  -  A1 ( 1 ,1) 
10     CONTINUE 
RETURN 
END 
C 
C 
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C    SUBROUTINE  STAT 
C 

SUBROUT INE  STAT ( N , R , BX , DX , BOX , DOX  ) 
INTEGER  J,R,N 

REAL  BX(N+1  ),DX(N+1  )  , BOX( R ) , DOX( R ) 
BOX(2)=  0.0 
DOX(2)=  0.0 
C 

C    COMPUTE  THE  ESTIMATES  OF  THE  SAMPLE  MEAN  AND 
C   VARIANCE 
C 

BOX( 1 )  =  BX(N+1 , J)/N 
DOX( 1 )  =  DX(N+1 , J)/N 
DO  260  1=1  ,N 

B0X(2)=  B0X(2)  +  (BX( I)-BOX( 1  )  )**2 
DOX(2)=  D0X(2)  +  (DX( I)-DOX( 1  )  )**2 
260        CONTINUE 

B0X(2)=  BOX(2)/(N*(N-1 )) 
D0X(2)=  D0X(2)/(N*(N-1  )  ) 
RETURN 
END 
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APPENDIX  D 

FORTRAN  PROGRAM  LISTING  OF  THE  BCD  SIMULATION  USING 

STRATIFIED  SAMPLING 


DIMENSION  BX( 101  )  ,DX( 101  )  ,SEED(2) ,B0X(2) ,DOX(2) 

INTEGER  I,N,BB,DD,R 

REAL*4  X, Y,Z, BX.DX.TXT, BOX, DOX 

DOUBLE  PRECISION  SEED 

DATA  SEED  / 1 234 . 0D0 , 5678901 23 . 0D0/ 

R=  2 

N=1  00 
C 

C   RECEIVE  INPUT  DATA  FROM  TERMINAL 
C 

WRITE( * ,3) 

3  FORMAT ( 1X, 'ENTER  THE  NUMBER  OF  BOMBERS') 
READ  ' (  12)  '  ,  BB 

WRITE( *  ,4) 

4  FORMAT ( 1X, 'ENTER  THE  RATE  WHICH  A  BOMBER  SHOOTS' 
&  'DOWN  A  DEFENDER') 

READ  ' (F5.3) ' ,  Y 
WRITE( *  ,5) 

5  FORMAT ( 1 X, 'ENTER  THE  NUMBER  OF  DEFENDERS') 
READ  '  ( 12)  '  ,  DD 

WRITE( * ,6) 

6  FORMAT ( 1X, 'ENTER  THE  RATE  WHICH  A  DEFENDER' 
&  'SHOOTS  DOWN  A  BOMBER') 

READ  ' (F5.3) ' ,  X 

WRITE( * ,7) 'ENTER  THE  RATE  WHICH   FREE  DEFENDERS' 
&  'FIND  FREE  BOMBERS' 

7  FORMAT (1X, A) 
READ  ' (F5.3) ' ,  Z 

WRITE( * ,8) 'ENTER  THE   TIME  DURATION  OF  THE' 
&  'BATTLE' 

8  FORMAT (1X, A) 

READ  ' (F5.3) ' ,  TXT 
C 

C    RUN  REPLICATIONS  OF  N  BATTLES  AND  OBTAIN  SUMMARY  OF 
C    N  BATTLES 
C 

CALL  BATTLE( N , R , SEED , X , Y , Z , BB , DD , TXT , BX , DX  ) 
C 

C   COMPUTE  STATISTICAL  OUTPUT  ANALYSIS  OBTAIN  PARAMETER 
C   ESTIMATES 
C 

CALL  STAT(N,R,BX,DX,BOX,DOX) 
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c 

C   FORMAT  AND  PRINT  OUTPUT  OF  PARAMETER  ESTIMATES 
C 

WRITE(3,279)  N 

279  FORMAT ( 1 5X , 'SAMPLE  SIZE ' , 16 , / 1 5X , 6 ( ' - ' ) , 
&  1X,4( '-' )) 

WRITE(3,280)  BB,DD 

280  &  FORMAT( //1X, 14, 2X, 'BOMBERS' ,6X, 'VERSUS' , 
&  6X, 14, 2X, 'DEFENDERS' , 

&  /1X, 13(  '-'  ), 18X, 14(  '-'  ) ) 
290    WRITE  (3,300)TXT 
300    FORMAT ( //18X,  'TIME'  ,F6. 1  ,/ ,47(  '-'  )  ) 

WRITE(3,310) 
3 1 0     FORMAT ( 1 9X , ' BOMBER ' , 2X , ' DEFENDER ' , 
&  /18X,7(  '  =  '  )  ,2X,8(  '  =  '  ) ) 

WRITE(3,320)  BOX( 1  )  ,DOX( 1  ) 
320    FORMAT( 1X, 'AVERAGE' ,7X,2F10.4) 

WRITE(3,330)  B0X(2) ,DOX(2) 
330    FORMAT( 1X, 'VARIANCE' ,6X,2F10.4) 
STOP 
END 
C 
C 

C    SUBROUTINE  BATTLE 
C 

SUBROUTINE 
&  BATTLE( N , R , SEED , X , Y , Z , BB , DD , TXT , BX , DX ) 
INTEGER  BB,DD,H,I,J,K,N,R,W 
REAL  GA(3, 1000) ,SOJ(3, 1000) ,BM(34,3) , 
&  DF(34,3) , BX(N+1 ) ,DX(N+1 ) , X , Y , Z , T , TXT , TIME , 
&  INF,NC,BW,DW,B,C,D,BM,DF 
DOUBLE  PRECISION  SEED(R) 
K=  1000 
W=  3 
C 

C  INITIALIZE  STATISTICAL  COUNTERS 
C 

BX(N+1 )=  0.0 
DX(N+1 )=  0.0 
C 

C  RUN  N  REPLICATIONS 
C 

DO  200  1=1  ,N 

CALL  SOJOUR(SEED( 1 ) ,SOJ,W,K) 
CALL  STATE(SEED(2) ,GA,W,K) 
DO  175  J=1  ,W 
C 

C    INITIALIZE  START-TO-BATTLE  VALUES 
C 

H  =  0 
T  =  0  .  0 
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C  =  0.0 

B=REAL(BB) 

D=REAL(DD) 

NC=  1  . 

0 

BW=  0. 

0 

DW=  0 

0 

C 

C   OBSERVE  AN  OCCURRENCE  OF  AN  EVENT 

C 

100              H=H+1 

C 

C  DETERMINE  NEXT  INTERIM  EVENT  AND  UPDATE  THE  STATE 

C       (B,C,D) 

C 

IF(GA( J 

,H)  .LE.  BW)  THEN 

B  =  B  + 

1  .0 

C  =  C  - 

1  .0 

D  =  D 

ELSE  IF 

( GA ( J , H )  . LE .  ( BW+NC ) )  THEN 

B  =  B  - 

1  .0 

C  =  C  + 

1  .0 

D  =  D  - 

1  .0 

ELSE 

B=B 

C  =  C  - 

1  .0 

D  =  D+ 

1  .0 

END  IF 

c 

C   COMPUTE  MEAN  TIME  IN  STATE  (B,C,D) 
C 

IF((B  .EQ.   0.0  .OR.  D  .EQ.  0 . 0  )  . AND .  C 
&  .EQ.  0.0)  THEN 

INF=  1000000.0 
ELSE 

INF=  1.0/(Z*B*D  +  (X+Y)*C) 
END  IF 
C 

C   COMPUTE  SOJOURN  TIME  OF  STATE  (B,C,D) 
C 

TIME=  -INF  *  ALOG(SOJ( J,H) ) 
C 

C   ADVANCE  SIMULATED  TIME  OF  THE  AIR  BATTLE 
C 

T=  T  +  TIME 
C 
C   COMPUTE  PROBABILITY  OF  NEXT  INTERIM  EVENTS  OCCURRING 


NC=  Z*B*D*INF 
BW=  Y*C*INF 
DW=  X*C*INF 
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c 

C   CHECK  FOR  OCCURRENCE  OF  TERMINATION  EVENT 
C 

IF  (T  .LT.  TXT(TI))  GOTO  100 
C 

C    RECORD  RESULTS  OF  BATTLE 
C 

BM( I, J)=  B  +  C 
DF( I, J)=  D  +  C 
175       CONTINUE 
C 

C    ACCUMULATE  SUMMARY  OF  N  BATTLE  OUTCOMES 
C 

BX(I)=  (BM(I,1)  +  BM(I,2)  +  BM(I,3)/3.0 
DX(I)=  (DF(I,1)  +  DF(I,2)  +  BM(I,3)/3.0 
BX(N+1  )=  BX(N+1  )  +  BX(I) 
DX(N+1  )=  DX(N+1  )  +  DX(I) 
200    CONTINUE 
RETURN 
END 
C 
C 

C    SUBROUTINE  SOJOUR 
C 

SUBROUTINE  SOJOUR( SEED2 , B2 , W , K ) 

INTEGER  I,W,K,J 

REAL  B2(W,K),A2 

DOUBLE  PRECISION  EFF , SEED2 

EFF=  2147483647. 0D0 

DO  10  1=1  ,K 

SEED2=DM0D( 1 6807.0D0  *  SEED2 , EFF ) 
A2=  SEED2/EFF 
DO  5  J-1  ,W 

B2(J,I)=  AM0D(A2  +  ((J-1)  *  1.0)/3.0,1 
5         CONTINUE 
10     CONTINUE 
RETURN 
END 
C 
C 

C    SUBROUTINE  STATE 
C 

SUBROUTINE  STATE( SEED1 , A1 , W , K ) 

INTEGER  I,K,W,J 

REAL  A1 (W,K) ,A2 

DOUBLE  PRECISION  EFF,SEED1 

EFF=  21 47483647. 0D0 

DO  10  1=1 ,K 

SEED1 =DMOD( 1 6807.0D0  *  SEED1,EFF) 
A2=  SEED 1 /EFF 
DO  5  J=1 ,W 
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A1(J,I)=  AM0D(A2  +  ((-J-1)  *  1 .0)/3.0, 1 .0) 
5         CONTINUE 
RETURN 
END 
C 
C 

C    SUBROUTINE  STAT 
C 

SUBROUTINE  STAT( N , R , BX , DX , BOX , DOX ) 
INTEGER  R,N 

REAL  BX(N+1  ),DX(N+1  )  , BOX ( R ) , DOX( R ) 
B0X(2)=  0.0 
D0X(2)=  0.0 
C 

C    COMPUTE  THE  ESTIMATES  OF  THE  SAMPLE  MEAN  AND 
C   VARIANCE 
C 

BOX( 1 )  =  BX(N+1 )/N 
DOX( 1 )  =  DX(N+1 )/N 
DO  260  1=1 ,N 

BOX(2)=  B0X(2)  +  (BX( I )-BOX( 1 ) )**2 
D0X(2)=  D0X(2)  +  (DX( I )-DOX( 1 ) )**2 
260        CONTINUE 

BOX(2)=  B0X(2)/(N«(N-1 )) 
D0X(2)=  DOX(2)/(N*(N-1 ) ) 
RETURN 
END 
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APPENDIX  E 
STATISTICAL  OUTPUT  DATA  FROM  THE  BCD  SIMULATIONS 


TABLE  E.1   OUTPUT  PARAMETERS  OF  THE  BCD  MODEL  ESTIMATED 
FROM  CRUDE  SIMULATION 


SCENARIO 

DEFENDER 

BOMBER 

SIMULATION 

DEFENDER 

BOMBER 

TIME 

MEAN 

VAR 

MEAN 

VAR 

1 

1  0 

10 

25 

7.6 

.01  91 

7.5 

.0250 

2 

10 

10 

75 

4.4 

.0387 

4.5 

.0358 

3 

1  0 

10 

125 

2.8 

.0375 

3.4 

.  041  7 

4 

1  0 

30 

25 

5.3 

.0328 

25.  1 

.0541 

5 

10 

30 

75 

1  .3 

.0157 

21  .3 

.  1  402 

6 

10 

30 

125 

.3 

.0036 

20.8 

.  1725 

7 

10 

50 

25 

4.2 

.0251 

44.7 

.0558 

8 

10 

50 

75 

.6 

.0064 

40.8 

.  1  654 

9 

10 

50 

125 

.  1 

.001  4 

39.7 

.  1874 

10 

30 

10 

25 

25.1 

.0458 

5.4 

.0350 

1  1 

30 

10 

75 

21  .  1 

.  1  407 

1  .4 

.0141 

12 

30 

10 

125 

19.8 

.2247 

.4 

.0045 

13 

30 

30 

25 

18.7 

.0933 

18.5 

.0842 

14 

30 

30 

75 

8.2 

.1618 

8.4 

.  1895 

15 

30 

30 

125 

5.7 

.  1984 

4.7 

.  1  341 

16 

30 

50 

25 

14.4 

.0628 

34.5 

.  1  1  04 

17 

30 

50 

75 

3.4 

.0362 

22.9 

.3479 

18 

30 

50 

125 

.9 

.0136 

20.  1 

.4881 

19 

50 

10 

25 

44.4 

.0651 

4.  1 

.0317 

20 

50 

10 

75 

40.7 

.  1427 

.7 

.0063 

21 

50 

10 

125 

41  .  1 

.  1666 

.  1 

.001  4 

22 

50 

30 

25 

34.3 

.  1687 

14.0 

.0826 

23 

50 

30 

75 

23.5 

.3339 

3.  1 

.0359 

24 

50 

30 

125 

21  .  1 

.5300 

.8 

.01  07 

25 

50 

50 

25 

27.9 

.  1590 

26.4 

.  1  501 

26 

50 

50 

75 

9.5 

.2717 

11.5 

.2865 

27 

50 

50 

125 

5.2 

.2126 

7.9 

.2986 
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TABLE  E.2   OUTPUT  PARAMETERS  OF  THE  BCD  MODEL  ESTIMATED 
FROM  THE  SIMULATION  USING  ANTITHETIC  VARIATES 


SCENARIO 

DEFENDER 

BOMBER 

SIMULATION 

DEFENDER  BOMBER 

TIME 

MEAN 

VAR 

MEAN 

VAR 

1 

10 

10 

25 

7.6 

.0067 

7.5 

.0105 

2 

10 

10 

75 

4.5 

.0048 

4.4 

.0060 

3 

10 

10 

125 

3.2 

.0079 

3.3 

.0057 

4 

10 

30 

25 

5.1 

.0129 

25  .  1 

.  0100 

5 

10 

30 

75 

1  .2 

.0071 

21  .3 

.0329 

6 

10 

30 

125 

.3 

.  0024 

20.2 

.  0446 

7 

10 

50 

25 

4.4 

.0123 

44.4 

.0146 

8 

10 

50 

75 

.6 

.0060 

40.9 

.0395 

9 

10 

50 

125 

.  1 

.0013 

40.  1 

.0683 

10 

30 

10 

25 

25.4 

.0116 

5.0 

.0116 

11 

30 

10 

75 

21.  1 

.  0271 

1.2 

.0056 

12 

30 

10 

125 

20.2 

.0558 

.3 

.0030 

13 

30 

30 

25 

18.5 

.0164 

18.4 

.0210 

14 

30 

30 

75 

8.5 

.0251 

8.4 

.0270 

15 

30 

30 

125 

5.  1 

.0307 

5  .  1 

.0352 

16 

30 

50 

25 

14.6 

.0223 

34.6 

.  0308 

17 

30 

50 

75 

3.0 

.0138 

23.6 

.0542 

18 

30 

50 

125 

1.0 

.0075 

20.6 

.0446 

19 

50 

10 

25 

44.0 

.0169 

4.3 

.0104 

20 

50 

10 

75 

40.2 

.0410 

.6 

.0059 

21 

50 

10 

125 

40.6 

.0451 

.  1 

.0015 

22 

50 

30 

25 

34.3 

.0279 

14.4 

.  0368 

23 

50 

30 

75 

23.1 

.  0375 

3.1 

.  0194 

24 

50 

30 

125 

20.5 

.0739 

.9 

.0077 

25 

50 

50 

25 

27.3 

.0266 

27.  1 

.  0357 

26 

50 

50 

75 

10.5 

.0269 

10.7 

.0309 

27 

50 

50 

125 

6.4 

.0458 

6.4 

.0472 
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TABLE  E.3 

EFFICIENCY 

OF  ANTITHETIC  VARIATES 

IN  THE  BCD 

MODEL 

SCENARIO 

DEFENDER 

BOMBER 

SIMULATION 

DEFENDER 

BOMBER 

TIME 

RE 

IP(*) 

RE 

IP(#) 

1 

10 

10 

25 

2.9 

65.0 

2.4 

58.  0 

2 

10 

10 

75 

8.  1 

87.6 

6.0 

83.2 

3 

10 

10 

125 

4.7 

78.9 

7.3 

86.3 

4 

10 

30 

25 

2.5 

60.7 

5  .4 

81  .  5 

5 

10 

30 

75 

2.2 

54  .8 

4  .3 

76.  5 

6 

10 

30 

125 

1.5 

33.3 

3.9 

74.1 

7 

10 

50 

25 

2.0 

51.0 

3.8 

73.8 

8 

10 

50 

75 

1.  1 

6.3 

4.2 

76.1 

9 

10 

50 

125 

1.1 

7.1 

2.7 

63.6 

10 

30 

10 

25 

3.9 

74.7 

3.0 

66.9 

11 

30 

10 

75 

5.2 

80.7 

2.5 

60.3 

12 

30 

10 

125 

4.0 

75.2 

1.5 

33.3 

13 

30 

30 

25 

5.7 

82.4 

4.0 

75.1 

14 

30 

30 

75 

6.4 

84.5 

7.0 

85.8 

15 

30 

30 

125 

6.5 

84.5 

3.8 

73.8 

16 

30 

50 

25 

2.8 

64.5 

3.6 

72.  1 

17 

30 

50 

75 

2.6 

61.9 

6.4 

84.4 

18 

30 

50 

125 

1.8 

44.9 

10.9 

90.9 

19 

50 

10 

25 

3.9 

74.0 

3.0 

67.2 

20 

50 

10 

75 

3.5 

71.2 

1.  1 

6.3 

21 

50 

10 

125 

3.7 

72.9 

.9 

7.1 

22 

50 

30 

25 

6.0 

83.4 

2.2 

55  .4 

23 

50 

30 

75 

8.9 

88.8 

1.9 

46.  0 

24 

50 

30 

125 

7.2 

86.1 

1.4 

28.0 

25 

50 

50 

25 

6.0 

83.3 

4.2 

76.2 

26 

50 

50 

75 

10.1 

90.1 

9.3 

89.2 

27 

50 

50 

125 

4.6 

78.5 

6.3 

84.2 
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TABLE  E.4   OUTPUT  PARAMETERS  OF  THE  BCD  MODEL  ESTIMATED 
FROM  SIMULATION  USING  STRATIFIED  SAMPLING 


SCENARIO 

DEFENDER 
MEAN    VAR 

BOMBE 
MEAN 

R 

SIMULATION 

DEFENDER 

BOMBER 

TIME 

VAR 

1 

10 

10 

25 

7.4 

.0119 

7.5 

.0114 

2 

10 

1  0 

75 

4.4 

.0248 

4.5 

.0237 

3 

10 

10 

125 

3.1 

.01  60 

3.2 

.01  97 

4 

10 

30 

25 

5.1 

.0124 

25.0 

.0252 

5 

10 

30 

75 

1  .2 

.0076 

21  .2 

.071  0 

6 

10 

30 

125 

.4 

.0038 

20.5 

.0801 

7 

10 

50 

25 

4.4 

.01  26 

44.0 

.0375 

8 

10 

50 

75 

.9 

.0083 

40.2 

.  1  021 

9 

10 

50 

125 

.  1 

.0006 

40  .0 

.0930 

1  0 

30 

10 

25 

25.  1 

.0233 

5.3 

.0124 

1  1 

30 

1  0 

75 

21  .  1 

.0726 

1  .3 

.0118 

1  2 

30 

10 

125 

20.3 

.0877 

.4 

.0027 

13 

30 

30 

25 

18.4 

.0283 

18.4 

.0193 

14 

30 

30 

75 

8.0 

.0439 

8.2 

.0700 

15 

30 

30 

125 

5.5 

.0505 

4.8 

.0426 

1  6 

30 

50 

25 

14.2 

.0219 

34.5 

.0798 

17 

30 

50 

75 

3.2 

.021  1 

23.3 

.  1236 

18 

30 

50 

125 

.8 

.0081 

20.  1 

.  1880 

19 

50 

10 

25 

44.6 

.0339 

4.3 

.0163 

20 

50 

10 

75 

41  .0 

.0872 

.7 

.0053 

21 

50 

10 

125 

40.4 

.  1  089 

.  1 

.0012 

22 

50 

30 

25 

34.6 

.0571 

14.3 

.0439 

23 

50 

30 

75 

22.6 

.  1  187 

3.3 

.0197 

24 

50 

30 

125 

20.9 

.2402 

1  .0 

.0086 

25 

50 

50 

25 

27.7 

.0630 

27.0 

.  1  084 

26 

50 

50 

75 

10.5 

.0898 

10.7 

.0881 

27 

50 

50 

125 

6.5 

.091  1 

6.4 

.  1  084 
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TABLE  E.5 

EFFICIENCY 

OF  STRATIFIED  SAMPLING 

IN  THE  BCD 

MODEL 

SCENARIO 

DEFENDER 
RE   IP(%) 

BOM 
RE 

BER 

SIMULATION 

DEFENDER 

BOMBER 

TIME 

IP(*) 

1 

10 

1  0 

25 

1  .6 

39.0 

2  .2 

54.4 

2 

10 

10 

75 

1  .6 

36.0 

1  .5 

33.8 

3 

10 

10 

125 

2.3 

57.3 

2.  1 

52.8 

4 

10 

30 

25 

2.6 

62.2 

2.  1 

53.4 

5 

10 

30 

75 

2.  1 

51  .6 

2.0 

49.4 

6 

10 

30 

125 

.9 

-5.5 

2.2 

53.6 

7 

10 

50 

25 

2.0 

49.8 

1  .5 

32  .8 

8 

10 

50 

75 

.8 

-29.7 

1  .6 

38.3 

9 

10 

50 

125 

2.3 

57.  1 

2.0 

50.4 

1  0 

30 

10 

25 

2.0 

49.  1 

2.8 

64.6 

1  1 

30 

10 

75 

1  .9 

48.4 

1  .2 

16.3 

1  2 

30 

10 

125 

2.6 

61  .0 

1  .7 

40.0 

13 

30 

30 

25 

3.3 

69.7 

4.4 

77.  1 

14 

30 

30 

75 

3.7 

72.9 

2.7 

63.  1 

15 

30 

30 

125 

3.9 

74.5 

3.1 

68.2 

16 

30 

50 

25 

2.9 

65.  1 

1  .4 

27.7 

17 

30 

50 

75 

1  .7 

41  .7 

2.8 

64.5 

18 

30 

50 

125 

1  .7 

40.7 

2.6 

61  .5 

19 

50 

1  0 

25 

1  .9 

47.9 

1  .9 

48.6 

20 

50 

10 

75 

1  .6 

39.0 

1  .2 

15.9 

21 

50 

10 

125 

1  .5 

34.6 

1  .2 

14.2 

22 

50 

30 

25 

3.0 

66.2 

1  .9 

46.9 

23 

50 

30 

75 

2.8 

64.5 

1  .8 

45.  1 

24 

50 

30 

125 

2.2 

54.7 

1  .2 

19.6 

25 

50 

50 

25 

2.5 

60.4 

1  .4 

27.8 

26 

50 

50 

75 

3.0 

66.9 

3.3 

69.2 

27 

50 

50 

125 

2.3 

57.  1 

2.8 

63.7 
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